Architecture and membrane interactions of the EGF receptor.

Dimerization-driven activation of the intracellular kinase domains of the epidermal growth factor receptor (EGFR) upon extracellular ligand binding is crucial to cellular pathways regulating proliferation, migration, and differentiation. Inactive EGFR can exist as both monomers and dimers, suggesting that the mechanism regulating EGFR activity may be subtle. The membrane itself may play a role but creates substantial difficulties for structural studies. Our molecular dynamics simulations of membrane-embedded EGFR suggest that, in ligand-bound dimers, the extracellular domains assume conformations favoring dimerization of the transmembrane helices near their N termini, dimerization of the juxtamembrane segments, and formation of asymmetric (active) kinase dimers. In ligand-free dimers, by holding apart the N termini of the transmembrane helices, the extracellular domains instead favor C-terminal dimerization of the transmembrane helices, juxtamembrane segment dissociation and membrane burial, and formation of symmetric (inactive) kinase dimers. Electrostatic interactions of EGFR's intracellular module with the membrane are critical in maintaining this coupling.

Oncogenic mutations counteract intrinsic disorder in the EGFR kinase and promote receptor dimerization.

The mutation and overexpression of the epidermal growth factor receptor (EGFR) are associated with the development of a variety of cancers, making this prototypical dimerization-activated receptor tyrosine kinase a prominent target of cancer drugs. Using long-timescale molecular dynamics simulations, we find that the N lobe dimerization interface of the wild-type EGFR kinase domain is intrinsically disordered and that it becomes ordered only upon dimerization. Our simulations suggest, moreover, that some cancer-linked mutations distal to the dimerization interface, particularly the widespread L834R mutation (also referred to as L858R), facilitate EGFR dimerization by suppressing this local disorder. Corroborating these findings, our biophysical experiments and kinase enzymatic assays indicate that the L834R mutation causes abnormally high activity primarily by promoting EGFR dimerization rather than by allowing activation without dimerization. We also find that phosphorylation of EGFR kinase domain at Tyr845 may suppress the intrinsic disorder, suggesting a molecular mechanism for autonomous EGFR signaling.

Accurate and efficient integration for molecular dynamics simulations at constant temperature and pressure.

In molecular dynamics simulations, control over temperature and pressure is typically achieved by augmenting the original system with additional dynamical variables to create a thermostat and a barostat, respectively. These variables generally evolve on timescales much longer than those of particle motion, but typical integrator implementations update the additional variables along with the particle positions and momenta at each time step. We present a framework that replaces the traditional integration procedure with separate barostat, thermostat, and Newtonian particle motion updates, allowing thermostat and barostat updates to be applied infrequently. Such infrequent updates provide a particularly substantial performance advantage for simulations parallelized across many computer processors, because thermostat and barostat updates typically require communication among all processors. Infrequent updates can also improve accuracy by alleviating certain sources of error associated with limited-precision arithmetic. In addition, separating the barostat, thermostat, and particle motion update steps reduces certain truncation errors, bringing the time-average pressure closer to its target value. Finally, this framework, which we have implemented on both general-purpose and special-purpose hardware, reduces software complexity and improves software modularity.

Principles of conduction and hydrophobic gating in K+ channels.

We present the first atomic-resolution observations of permeation and gating in a K(+) channel, based on molecular dynamics simulations of the Kv1.2 pore domain. Analysis of hundreds of simulated permeation events revealed a detailed conduction mechanism, resembling the Hodgkin-Keynes "knock-on" model, in which translocation of two selectivity filter-bound ions is driven by a third ion; formation of this knock-on intermediate is rate determining. In addition, at reverse or zero voltages, we observed pore closure by a novel "hydrophobic gating" mechanism: A dewetting transition of the hydrophobic pore cavity-fastest when K(+) was not bound in selectivity filter sites nearest the cavity-caused the open, conducting pore to collapse into a closed, nonconducting conformation. Such pore closure corroborates the idea that voltage sensors can act to prevent pore collapse into the intrinsically more stable, closed conformation, and it further suggests that molecular-scale dewetting facilitates a specific biological function: K(+) channel gating. Existing experimental data support our hypothesis that hydrophobic gating may be a fundamental principle underlying the gating of voltage-sensitive K(+) channels. We suggest that hydrophobic gating explains, in part, why diverse ion channels conserve hydrophobic pore cavities, and we speculate that modulation of cavity hydration could enable structural determination of both open and closed channels.

Refinement of protein structure homology models via long, all-atom molecular dynamics simulations.

Accurate computational prediction of protein structure represents a longstanding challenge in molecular biology and structure-based drug design. Although homology modeling techniques are widely used to produce low-resolution models, refining these models to high resolution has proven difficult. With long enough simulations and sufficiently accurate force fields, molecular dynamics (MD) simulations should in principle allow such refinement, but efforts to refine homology models using MD have for the most part yielded disappointing results. It has thus far been unclear whether MD-based refinement is limited primarily by accessible simulation timescales, force field accuracy, or both. Here, we examine MD as a technique for homology model refinement using all-atom simulations, each at least 100 µs long-more than 100 times longer than previous refinement simulations-and a physics-based force field that was recently shown to successfully fold a structurally diverse set of fast-folding proteins. In MD simulations of 24 proteins chosen from the refinement category of recent Critical Assessment of Structure Prediction (CASP) experiments, we find that in most cases, simulations initiated from homology models drift away from the native structure. Comparison with simulations initiated from the native structure suggests that force field accuracy is the primary factor limiting MD-based refinement. This problem can be mitigated to some extent by restricting sampling to the neighborhood of the initial model, leading to structural improvement that, while limited, is roughly comparable to the leading alternative methods.

A conserved protonation-dependent switch controls drug binding in the Abl kinase.

In many protein kinases, a characteristic conformational change (the "DFG flip") connects catalytically active and inactive conformations. Many kinase inhibitors--including the cancer drug imatinib--selectively target a specific DFG conformation, but the function and mechanism of the flip remain unclear. Using long molecular dynamics simulations of the Abl kinase, we visualized the DFG flip in atomic-level detail and formulated an energetic model predicting that protonation of the DFG aspartate controls the flip. Consistent with our model's predictions, we demonstrated experimentally that the kinetics of imatinib binding to Abl kinase have a pH dependence that disappears when the DFG aspartate is mutated. Our model suggests a possible explanation for the high degree of conservation of the DFG motif: that the flip, modulated by electrostatic changes inherent to the catalytic cycle, allows the kinase to access flexible conformations facilitating nucleotide binding and release.

How does a drug molecule find its target binding site?

Although the thermodynamic principles that control the binding of drug molecules to their protein targets are well understood, detailed experimental characterization of the process by which such binding occurs has proven challenging. We conducted relatively long, unguided molecular dynamics simulations in which a ligand (the cancer drug dasatinib or the kinase inhibitor PP1) was initially placed at a random location within a box that also contained a protein (Src kinase) to which that ligand was known to bind. In several of these simulations, the ligand correctly identified its target binding site, forming a complex virtually identical to the crystallographically determined bound structure. The simulated trajectories provide a continuous, atomic-level view of the entire binding process, revealing persistent and noteworthy intermediate conformations and shedding light on the role of water molecules. The technique we employed, which does not assume any prior knowledge of the binding site's location, may prove particularly useful in the development of allosteric inhibitors that target previously undiscovered binding sites.

Dynamic control of slow water transport by aquaporin 0: implications for hydration and junction stability in the eye lens.

Aquaporin 0 (AQP0), the most abundant membrane protein in mammalian lens fiber cells, not only serves as the primary water channel in this tissue but also appears to mediate the formation of thin junctions between fiber cells. AQP0 is remarkably less water permeable than other aquaporins, but the structural basis and biological significance of this low permeability remain uncertain, as does the permeability of the protein in a reported junctional form. To address these issues, we performed molecular dynamics (MD) simulations of water transport through membrane-embedded AQP0 in both its (octameric) junctional and (tetrameric) nonjunctional forms. From our simulations, we measured an osmotic permeability for the nonjunctional form that agrees with experiment and found that the distinct dynamics of the conserved, lumen-protruding side chains of Tyr-23 and Tyr-149 modulate water passage, accounting for the slow permeation. The junctional and nonjunctional forms conducted water equivalently, in contrast to a previous suggestion based on static crystal structures that water conduction is lost on junction formation. Our analysis suggests that the low water permeability of AQP0 may help maintain the mechanical stability of the junction. We hypothesize that the structural features leading to low permeability may have evolved in part to allow AQP0 to form junctions that both conduct water and contribute to the organizational structure of the fiber cell tissue and microcirculation within it, as required to maintain transparency of the lens.

Microsecond molecular dynamics simulation shows effect of slow loop dynamics on backbone amide order parameters of proteins.

A molecular-level understanding of the function of a protein requires knowledge of both its structural and dynamic properties. NMR spectroscopy allows the measurement of generalized order parameters that provide an atomistic description of picosecond and nanosecond fluctuations in protein structure. Molecular dynamics (MD) simulation provides a complementary approach to the study of protein dynamics on similar time scales. Comparisons between NMR spectroscopy and MD simulations can be used to interpret experimental results and to improve the quality of simulation-related force fields and integration methods. However, apparent systematic discrepancies between order parameters extracted from simulations and experiments are common, particularly for elements of noncanonical secondary structure. In this paper, results from a 1.2 micros explicit solvent MD simulation of the protein ubiquitin are compared with previously determined backbone order parameters derived from NMR relaxation experiments [Tjandra, N.; Feller, S. E.; Pastor, R. W.; Bax, A. J. Am. Chem. Soc. 1995, 117, 12562-12566]. The simulation reveals fluctuations in three loop regions that occur on time scales comparable to or longer than that of the overall rotational diffusion of ubiquitin and whose effects would not be apparent in experimentally derived order parameters. A coupled analysis of internal and overall motion yields simulated order parameters substantially closer to the experimentally determined values than is the case for a conventional analysis of internal motion alone. Improved agreement between simulation and experiment also is encouraging from the viewpoint of assessing the accuracy of long MD simulations.

Assessment of the utility of contact-based restraints in accelerating the prediction of protein structure using molecular dynamics simulations.

Molecular dynamics (MD) simulation is a well-established tool for the computational study of protein structure and dynamics, but its application to the important problem of protein structure prediction remains challenging, in part because extremely long timescales can be required to reach the native structure. Here, we examine the extent to which the use of low-resolution information in the form of residue-residue contacts, which can often be inferred from bioinformatics or experimental studies, can accelerate the determination of protein structure in simulation. We incorporated sets of 62, 31, or 15 contact-based restraints in MD simulations of ubiquitin, a benchmark system known to fold to the native state on the millisecond timescale in unrestrained simulations. One-third of the restrained simulations folded to the native state within a few tens of microseconds-a speedup of over an order of magnitude compared with unrestrained simulations and a demonstration of the potential for limited amounts of structural information to accelerate structure determination. Almost all of the remaining ubiquitin simulations reached near-native conformations within a few tens of microseconds, but remained trapped there, apparently due to the restraints. We discuss potential methodological improvements that would facilitate escape from these near-native traps and allow more simulations to quickly reach the native state. Finally, using a target from the Critical Assessment of protein Structure Prediction (CASP) experiment, we show that distance restraints can improve simulation accuracy: In our simulations, restraints stabilized the native state of the protein, enabling a reasonable structural model to be inferred.

Atomic-level simulation of current-voltage relationships in single-file ion channels.

The difficulty in characterizing ion conduction through membrane channels at the level of individual permeation events has made it challenging to elucidate the mechanistic principles underpinning this fundamental physiological process. Using long, all-atom simulations enabled by special-purpose hardware, we studied K(+) permeation across the KV1.2/2.1 voltage-gated potassium channel. At experimentally accessible voltages, which include the physiological range, the simulated permeation rate was substantially lower than the experimentally observed rate. The current-voltage relationship was also nonlinear but became linear at much higher voltages. We observed permeation consistent with a "knock-on" mechanism at all voltages. At high voltages, the permeation rate was in accordance with our previously reported KV1.2 pore-only simulations, after the simulated voltages from the previous study were recalculated using the correct method, new insight into which is provided here. Including the voltage-sensing domains in the simulated channel brought the linear current-voltage regime closer to the experimentally accessible voltages. The simulated permeation rate, however, still underestimated the experimental rate, because formation of the knock-on intermediate occurred too infrequently. Reducing the interaction strength between the ion and the selectivity filter did not increase conductance. In complementary simulations of gramicidin A, similar changes in interaction strength did increase the observed permeation rate. Permeation nevertheless remained substantially below the experimental value, largely because of infrequent ion recruitment into the pore lumen. Despite the need to apply large voltages to simulate the permeation process, the apparent voltage insensitivity of the permeation mechanism suggests that the direct simulation of permeation at the single-ion level can provide fundamental physiological insight into ion channel function. Notably, our simulations suggest that the knock-on permeation mechanisms in KV1.2 and KcsA may be different.

Rotational relaxation in ortho-terphenyl: using atomistic simulations to bridge theory and experiment.

Understanding the nature of the glass transition--the dramatic slowing of dynamics and eventual emergence of a disordered solid from a cooling liquid--is a fundamental challenge in physical science. A central characteristic of glass-forming liquids is a non-exponential main relaxation process. The extent of deviation from exponential relaxation typically becomes more pronounced on cooling. Theories that predict a growth of spatially heterogeneous dynamics as temperature is lowered can explain these observations. In apparent contradiction to these theories, however, some experiments suggest that certain substances--notably including the intensely studied molecular glass-former ortho-terphenyl (OTP)--have a main relaxation process whose shape is essentially temperature independent, even though other observables predicted to be correlated with the degree of dynamical heterogeneity are temperature dependent. Here we report the first simulations based on an atomistic model of OTP that reach equilibrium at temperatures well into the supercooled regime. We first show that the results of these simulations are in reasonable quantitative agreement with experimental data for several basic properties over a wide range of temperatures. We then focus on rotational relaxation, finding nearly exponential behavior at high temperatures with clearly increasing deviations as temperature is lowered. The much weaker temperature dependence observed in light-scattering experiments also emerges from the same simulation data when we calculate correlation functions similar to those probed experimentally; this highlights the diversity of temperature dependencies that can be obtained with different probes. Further analysis suggests that the temperature insensitivity observed in the light-scattering experiments stems from the dependence of these measurements on internal as well as rotational molecular motion. Within the temperature range of our OTP simulations, our results strongly suggest that this archetypal glass-former behaves as anticipated by theories of the glass transition that predict increasing non-exponentiality with cooling, and our simulations thus strengthen the evidence supporting such theories.

Systematic validation of protein force fields against experimental data.

Molecular dynamics simulations provide a vehicle for capturing the structures, motions, and interactions of biological macromolecules in full atomic detail. The accuracy of such simulations, however, is critically dependent on the force field--the mathematical model used to approximate the atomic-level forces acting on the simulated molecular system. Here we present a systematic and extensive evaluation of eight different protein force fields based on comparisons of experimental data with molecular dynamics simulations that reach a previously inaccessible timescale. First, through extensive comparisons with experimental NMR data, we examined the force fields' abilities to describe the structure and fluctuations of folded proteins. Second, we quantified potential biases towards different secondary structure types by comparing experimental and simulation data for small peptides that preferentially populate either helical or sheet-like structures. Third, we tested the force fields' abilities to fold two small proteins--one α-helical, the other with ß-sheet structure. The results suggest that force fields have improved over time, and that the most recent versions, while not perfect, provide an accurate description of many structural and dynamical properties of proteins.

Equipartition and the Calculation of Temperature in Biomolecular Simulations.

Since the behavior of biomolecules can be sensitive to temperature, the ability to accurately calculate and control the temperature in molecular dynamics (MD) simulations is important. Standard analysis of equilibrium MD simulations-even constant-energy simulations with negligible long-term energy drift-often yields different calculated temperatures for different motions, however, in apparent violation of the statistical mechanical principle of equipartition of energy. Although such analysis provides a valuable warning that other simulation artifacts may exist, it leaves the actual value of the temperature uncertain. We observe that Tolman's generalized equipartition theorem should hold for long stable simulations performed using velocity-Verlet or other symplectic integrators, because the simulated trajectory is thought to sample almost exactly from a continuous trajectory generated by a shadow Hamiltonian. From this we conclude that all motions should share a single simulation temperature, and we provide a new temperature estimator that we test numerically in simulations of a diatomic fluid and of a solvated protein. Apparent temperature variations between different motions observed using standard estimators do indeed disappear when using the new estimator. We use our estimator to better understand how thermostats and barostats can exacerbate integration errors. In particular, we find that with large (albeit widely used) time steps, the common practice of using two thermostats to remedy so-called hot solvent-cold solute problems can have the counterintuitive effect of causing temperature imbalances. Our results, moreover, highlight the utility of multiple-time step integrators for accurate and efficient simulation.

Atomic-level characterization of the structural dynamics of proteins.

Molecular dynamics (MD) simulations are widely used to study protein motions at an atomic level of detail, but they have been limited to time scales shorter than those of many biologically critical conformational changes. We examined two fundamental processes in protein dynamics--protein folding and conformational change within the folded state--by means of extremely long all-atom MD simulations conducted on a special-purpose machine. Equilibrium simulations of a WW protein domain captured multiple folding and unfolding events that consistently follow a well-defined folding pathway; separate simulations of the protein's constituent substructures shed light on possible determinants of this pathway. A 1-millisecond simulation of the folded protein BPTI reveals a small number of structurally distinct conformational states whose reversible interconversion is slower than local relaxations within those states by a factor of more than 1000.

Minimizing thermodynamic length to select intermediate states for free-energy calculations and replica-exchange simulations.

In computational thermodynamics, a sequence of intermediate states is commonly introduced to connect two equilibrium states. We consider two cases where the choice of intermediate states is particularly important: minimizing statistical error in free-energy difference calculations and maximizing average acceptance probabilities in replica-exchange simulations. We derive bounds for these quantities in terms of the thermodynamic distance between the intermediates, and show that in both cases the intermediates should be chosen as equidistant points along a geodesic connecting the end states.

Automated Event Detection and Activity Monitoring in Long Molecular Dynamics Simulations.

Events of scientific interest in molecular dynamics (MD) simulations, including conformational changes, folding transitions, and translocations of ligands and reaction products, often correspond to high-level structural rearrangements that alter contacts between molecules or among different parts of a molecule. Due to advances in computer architecture and software, MD trajectories representing such structure-changing events have become easier to generate, but the length of these trajectories poses a challenge to scientific interpretation and analysis. In this paper, we present automated methods for the detection of potentially important structure-changing events in long MD trajectories. In contrast with traditional tools for the analysis of such trajectories, our methods provide a detailed report of broken and formed contacts that aids in the identification of specific time-dependent side-chain interactions. Our approach employs a coarse-grained representation of amino acid side chains, a contact metric based on higher order generalizations of Delaunay tetrahedralization, techniques for detecting significant shifts in the resulting contact time series, and a new kernel-based measure of contact alteration activity. The analysis methods we describe are incorporated in a newly developed package, called TimeScapes, which is freely available and compatible with trajectories generated by a variety of popular MD programs. Tests based on actual microsecond time scale simulations demonstrate that the package can be used to efficiently detect and characterize important conformational changes in realistic protein systems.

Mechanism of Na+/H+ antiporting.

Na+/H+ antiporters are central to cellular salt and pH homeostasis. The structure of Escherichia coli NhaA was recently determined, but its mechanisms of transport and pH regulation remain elusive. We performed molecular dynamics simulations of NhaA that, with existing experimental data, enabled us to propose an atomically detailed model of antiporter function. Three conserved aspartates are key to our proposed mechanism: Asp164 (D164) is the Na+-binding site, D163 controls the alternating accessibility of this binding site to the cytoplasm or periplasm, and D133 is crucial for pH regulation. Consistent with experimental stoichiometry, two protons are required to transport a single Na+ ion: D163 protonates to reveal the Na+-binding site to the periplasm, and subsequent protonation of D164 releases Na+. Additional mutagenesis experiments further validated the model.

Simulation studies of the fidelity of biomolecular structure ensemble recreation.

We examine the ability of Bayesian methods to recreate structural ensembles for partially folded molecules from averaged data. Specifically we test the ability of various algorithms to recreate different transition state ensembles for folding proteins using a multiple replica simulation algorithm using input from "gold standard" reference ensembles that were first generated with a Go-like Hamiltonian having nonpairwise additive terms. A set of low resolution data, which function as the "experimental" phi values, were first constructed from this reference ensemble. The resulting phi values were then treated as one would treat laboratory experimental data and were used as input in the replica reconstruction algorithm. The resulting ensembles of structures obtained by the replica algorithm were compared to the gold standard reference ensemble, from which those "data" were, in fact, obtained. It is found that for a unimodal transition state ensemble with a low barrier, the multiple replica algorithm does recreate the reference ensemble fairly successfully when no experimental error is assumed. The Kolmogorov-Smirnov test as well as principal component analysis show that the overlap of the recovered and reference ensembles is significantly enhanced when multiple replicas are used. Reduction of the multiple replica ensembles by clustering successfully yields subensembles with close similarity to the reference ensembles. On the other hand, for a high barrier transition state with two distinct transition state ensembles, the single replica algorithm only samples a few structures of one of the reference ensemble basins. This is due to the fact that the phi values are intrinsically ensemble averaged quantities. The replica algorithm with multiple copies does sample both reference ensemble basins. In contrast to the single replica case, the multiple replicas are constrained to reproduce the average phi values, but allow fluctuations in phi for each individual copy. These fluctuations facilitate a more faithful sampling of the reference ensemble basins. Finally, we test how robustly the reconstruction algorithm can function by introducing errors in phi comparable in magnitude to those suggested by some authors. In this circumstance we observe that the chances of ensemble recovery with the replica algorithm are poor using a single replica, but are improved when multiple copies are used. A multimodal transition state ensemble, however, turns out to be more sensitive to large errors in phi (if appropriately gauged) and attempts at successful recreation of the reference ensemble with simple replica algorithms can fall short.

Protein Structure Prediction: The Next Generation.

Over the last 10-15 years a general understanding of the chemical reaction of protein folding has emerged from statistical mechanics. The lessons learned from protein folding kinetics based on energy landscape ideas have benefited protein structure prediction, in particular the development of coarse grained models. We survey results from blind structure prediction. We explore how second generation prediction energy functions can be developed by introducing information from an ensemble of previously simulated structures. This procedure relies on the assumption of a funneled energy landscape keeping with the principle of minimal frustration. First generation simulated structures provide an improved input for associative memory energy functions in comparison to the experimental protein structures chosen on the basis of sequence alignment.