Atomistic view of the conformational activation of Src kinase using the string method with swarms-of-trajectories

CHARMM at 45: Enhancements in Accessibility, Functionality, and Speed

Since its inception nearly a half century ago, CHARMM has been playing a central role in computational biochemistry and biophysics. Commensurate with the developments in experimental research and advances in computer hardware, the range of methods and applicability of CHARMM have also grown. This review summarizes major developments that occurred after 2009 when the last review of CHARMM was published. They include the following: new faster simulation engines, accessible user interfaces for convenient workflows, and a vast array of simulation and analysis methods that encompass quantum mechanical, atomistic, and coarse-grained levels, as well as extensive coverage of force fields. In addition to providing the current snapshot of the CHARMM development, this review may serve as a starting point for exploring relevant theories and computational methods for tackling contemporary and emerging problems in biomolecular systems. CHARMM is freely available for academic and nonprofit research at https://academiccharmm.org/program.

Computational Analysis of Activation of Dimerized Epidermal Growth Factor Receptor Kinase Using the String Method and Markov State Model

Epidermal growth factor receptor (EGFR) activation is accompanied by dimerization. During the activation of the intracellular kinase domain, two EGFR kinases form an asymmetric dimer, and one side of the dimer (receiver) is activated. Using the string method and Markov state model (MSM), we performed a computational analysis of the structural changes in the activation of the EGFR dimer in this study. The string method reveals the minimum free-energy pathway (MFEP) from the inactive to active structure. The MSM was constructed from numerous trajectories of molecular dynamics simulations around the MFEP, which revealed the free-energy map of structural changes. In the activation of the receiver kinase, the unfolding of the activation loop (A-loop) is followed by the rearrangement of the C-helix, as observed in other kinases. However, unlike other kinases, the free-energy map of EGFR at the asymmetric dimer showed that the active state yielded the highest stability and revealed how interactions at the dimer interface induced receiver activation. As the H-helix of the activator approaches the C-helix of the receiver during activation, the A-loop unfolds. Subsequently, L782 of the receiver enters the pocket between the G- and H-helices of the activator, leading to a rearrangement of the hydrophobic residues around L782 of the receiver, which constitutes a structural rearrangement of the C-helix of the receiver from an outward to an inner position. The MSM analysis revealed long-time scale trajectories via kinetic Monte Carlo.

SHP-1 tyrosine phosphatase binding to c-Src kinase phosphor-dependent conformations: A comparative structural framework

SHP-1 is a cytosolic tyrosine phosphatase that is primarily expressed in hematopoietic cells. It acts as a negative regulator of numerous signaling pathways and controls multiple cellular functions involved in cancer pathogenesis. This study describes the binding preferences of SHP-1 (pY536) to c-Srcopen (pY416) and c-Srcclose (pY527) through in silico approaches. Molecular dynamics simulation analysis revealed more conformational changes in c-Srcclose upon binding to SHP-1, as compared to its active/open conformation that is stabilized by the cooperative binding of the C-SH2 domain and C-terminal tail of SHP-1 to c-Src SH2 and KD. In contrast, c-Srcclose and SHP-1 interaction is mediated by PTP domain-specific WPD-loop (WPDXGXP) and Q-loop (QTXXQYXF) binding to c-Srcclose C-terminal tail residues. The dynamic correlation analysis demonstrated a positive correlation for SHP-1 PTP with KD, SH3, and the C-terminal tail of c-Srcclose. In the case of the c-Srcopen-SHP-1 complex, SH3 and SH2 domains of c-Srcopen were correlated to C-SH2 and the C-terminal tail of SHP-1. Our findings reveal that SHP1-dependent c-Src activation through dephosphorylation relies on the conformational shift in the inhibitory C-terminal tail that may ease the recruitment of the N-SH2 domain to phosphotyrosine residue, resulting in the relieving of the PTP domain. Collectively, this study delineates the intermolecular interaction paradigm and underlying conformational readjustments in SHP-1 due to binding with the c-Src active and inactive state. This study will largely help in devising novel therapeutic strategies for targeting cancer development.

p38α Kinase Auto-Activation through Its Conformational Transition Induced by Tyr323 Phosphorylation

p38α is a key serine/threonine kinase that can enable atypical auto-activation through Zap70 phosphorylation and initiate T cell receptor signaling. The auto-activation plays an important role in autoimmune diseases. Although the classical activation mechanism of p38α has been studied in-depth, the atypical activation mechanism of Y323 phosphorylation-induced p38α auto-activation remains largely unexplained, especially the regulatory effects of phosphorylation on different sites (Y323 vs T180). From the X-ray experimental data, we identified the inactive and active states of p38α using principal component analysis. To understand the auto-activation process and the internal driving mechanism, a computational paradigm that couples the targeted molecular dynamics simulations, the String Method, and the umbrella sampling strategy were employed to generate the conformational landscape of p38α, including p38α T180-Y323, p38α T180-pY323, and p38α pT180-pY323 systems (pT180/pY323: phosphorylated T180/Y323). We explored that pY323 could change the conformational distribution and promote the conformational transition of p38α from the inactive state to the active state. Auto-activation of p38α is regulated by pY323 through destabilization of the hydrophobic core structure and aided by R173. This study will further explain the conformational transition of p38α induced by Y323 phosphorylation and provide insights into the universal molecular auto-activation mechanism of the p38 subfamily at the atomic level.

Identifying a Feasible Transition Pathway between Two Conformational States for a Protein

Proteins usually need to transit between different conformational states to fulfill their biological functions. In the mechanistic study of such transition processes by molecular dynamics simulations, identification of the minimum free energy path (MFEP) can substantially reduce the sampling space, thus enabling rigorous thermodynamic evaluation of the process. Conventionally, the MFEP is derived by iterative local optimization from an initial path, which is typically generated by simple brute force techniques like the targeted molecular dynamics (tMD). Therefore, the quality of the initial path determines the successfulness of MFEP estimation. In this work, we propose a method to improve derivation of the initial path. Through iterative relaxation-biasing simulations in a bidirectional manner, this method can construct a feasible transition pathway connecting two known states for a protein. Evaluation on small, fast-folding proteins against long equilibrium trajectories supports the good sampling efficiency of our method. When applied to larger proteins including the catalytic domain of human c-Src kinase as well as the converter domain of myosin VI, the paths generated by our method deviate significantly from those computed with the generic tMD approach. More importantly, free energy profiles and intermediate states obtained from our paths exhibit remarkable improvements over those from tMD paths with respect to both physical rationality and consistency with a priori knowledge.

A Companion Guide to the String Method with Swarms of Trajectories: Characterization, Performance, and Pitfalls

The string method with swarms of trajectories (SMwST) is an algorithm that identifies a physically meaningful transition pathway─a one-dimensional curve, embedded within a high-dimensional space of selected collective variables. The SMwST algorithm leans on a series of short, unbiased molecular dynamics simulations spawned at different locations of the discretized path, from whence an average dynamic drift is determined to evolve the string toward an optimal pathway. However conceptually simple in both its theoretical formulation and practical implementation, the SMwST algorithm is computationally intensive and requires a careful choice of parameters for optimal cost-effectiveness in applications to challenging problems in chemistry and biology. In this contribution, the SMwST algorithm is presented in a self-contained manner, discussing with a critical eye its theoretical underpinnings, applicability, inherent limitations, and use in the context of path-following free-energy calculations and their possible extension to kinetics modeling. Through multiple simulations of a prototypical polypeptide, combining the search of the transition pathway and the computation of the potential of mean force along it, several practical aspects of the methodology are examined with the objective of optimizing the computational effort, yet without sacrificing accuracy. In light of the results reported here, we propose some general guidelines aimed at improving the efficiency and reliability of the computed pathways and free-energy profiles underlying the conformational transitions at hand.

Exploring the Regulatory Function of the N-terminal Domain of SARS-CoV-2 Spike Protein through Molecular Dynamics Simulation

SARS-CoV-2 is what has caused the COVID-19 pandemic. Early viral infection is mediated by the SARS-CoV-2 homo-trimeric Spike (S) protein with its receptor binding domains (RBDs) in the receptor-accessible state. Molecular dynamics simulation on the S protein with a focus on the function of its N-terminal domains (NTDs) is performed. The study reveals that the NTD acts as a "wedge" and plays a crucial regulatory role in the conformational changes of the S protein. The complete RBD structural transition is allowed only when the neighboring NTD that typically prohibits the RBD's movements as a wedge detaches and swings away. Based on this NTD "wedge" model, it is proposed that the NTD-RBD interface should be a potential drug target.

TYROSINE KINASES: COMPLEX MOLECULAR SYSTEMS CHALLENGING COMPUTATIONAL METHODOLOGIES

Classical molecular dynamics (MD) simulations based on atomic models play an increasingly important role in a wide range of applications in physics, biology, and chemistry. Nonetheless, generating genuine knowledge about biological systems using MD simulations remains challenging. Protein tyrosine kinases are important cellular signaling enzymes that regulate cell growth, proliferation, metabolism, differentiation, and migration. Due to the large conformational changes and long timescales involved in their function, these kinases present particularly challenging problems to modern computational and theoretical frameworks aimed at elucidating the dynamics of complex biomolecular systems. Markov state models have achieved limited success in tackling the broader conformational ensemble and biased methods are often employed to examine specific long timescale events. Recent advances in machine learning continue to push the limitations of current methodologies and provide notable improvements when integrated with the existing frameworks. A broad perspective is drawn from a critical review of recent studies.

String Method with Swarms-of-Trajectories, Mean Drifts, Lag Time, and Committor

The kinetics of a dynamical system comprising two metastable states is formulated in terms of a finite-time propagator in phase space (position and velocity) adapted to the underdamped Langevin equation. Dimensionality reduction to a subspace of collective variables yields familiar expressions for the propagator, committor, and steady-state flux. A quadratic expression for the steady-state flux between the two metastable states can serve as a robust variational principle to determine an optimal approximate committor expressed in terms of a set of collective variables. The theoretical formulation is exploited to clarify the foundation of the string method with swarms-of-trajectories, which relies on the mean drift of short trajectories to determine the optimal transition pathway. It is argued that the conditions for Markovity within a subspace of collective variables may not be satisfied with an arbitrary short time-step and that proper kinetic behaviors appear only when considering the effective propagator for longer lag times. The effective propagator with finite lag time is amenable to an eigenvalue-eigenvector spectral analysis, as elaborated previously in the context of position-based Markov models. The time-correlation functions calculated by swarms-of-trajectories along the string pathway constitutes a natural extension of these developments. The present formulation provides a powerful theoretical framework to characterize the optimal pathway between two metastable states of a system.

Assessing the Performance of Traveling-salesman based Automated Path Searching (TAPS) on Complex Biomolecular Systems

Though crucial for understanding the function of large biomolecular systems, locating the minimum free energy paths (MFEPs) between their key conformational states is far from trivial due to their high-dimensional nature. Most existing path-searching methods require a static collective variable space as input, encoding intuition or prior knowledge of the transition mechanism. Such information is, however, hardly available a priori and expensive to validate. To alleviate this issue, we have previously introduced a Traveling-salesman based Automated Path Searching method (TAPS) and demonstrated its efficiency on simple peptide systems. Having implemented a parallel version of this method, here we assess the performance of TAPS on three realistic systems (tens to hundreds of residues) in explicit solvents. We show that TAPS successfully located the MFEP for the ground/excited state transition of the T4 lysozyme L99A variant, consistent with previous findings. TAPS also helped identifying the important role of the two polar contacts in directing the loop-in/loop-out transition of the mitogen-activated protein kinase kinase (MEK1), which explained previous mutant experiments. Remarkably, at a minimal cost of 126 ns sampling, TAPS revealed that the Ltn40/Ltn10 transition of lymphotactin needs no complete unfolding/refolding of its β-sheets and that five polar contacts are sufficient to stabilize the various partially unfolded intermediates along the MFEP. These results present TAPS as a general and promising tool for studying the functional dynamics of complex biomolecular systems.

The Perturbed Free-Energy Landscape: Linking Ligand Binding to Biomolecular Folding

The effects of ligand binding on biomolecular conformation are crucial in drug design, enzyme mechanisms, the regulation of gene expression, and other biological processes. Descriptive models such as "lock and key", "induced fit", and "conformation selection" are common ways to interpret such interactions. Another historical model, linked equilibria, proposes that the free-energy landscape (FEL) is perturbed by the addition of ligand binding energy for the bound population of biomolecules. This principle leads to a unified, quantitative theory of ligand-induced conformation change, building upon the FEL concept. We call the map of binding free energy over biomolecular conformational space the "binding affinity landscape" (BAL). The perturbed FEL predicts/explains ligand-induced conformational changes conforming to all common descriptive models. We review recent experimental and computational studies that exemplify the perturbed FEL, with emphasis on RNA. This way of understanding ligand-induced conformation dynamics motivates new experimental and theoretical approaches to ligand design, structural biology and systems biology.

Increasing the Sampling Efficiency of Protein Conformational Change by Combining a Modified Replica Exchange Molecular Dynamics and Normal Mode Analysis

Understanding conformational change at an atomic level is significant when determining a protein functional mechanism. Replica exchange molecular dynamics (REMD) is a widely used enhanced sampling method to explore protein conformational space. However, REMD with an explicit solvent model requires huge computational resources, immensely limiting its application. In this study, a variation of parallel tempering metadynamics (PTMetaD) with the omission of solvent-solvent interactions in exchange attempts and the use of low-frequency modes calculated by normal-mode analysis (NMA) as collective variables (CVs), namely ossPTMetaD, is proposed with the aim to accelerate MD simulations simultaneously in temperature and geometrical spaces. For testing the performance of ossPTMetaD, five protein systems with diverse biological functions and motion patterns were selected, including large-scale domain motion (AdK), flap movement (HIV-1 protease and BACE1), and DFG-motif flip in kinases (p38α and c-Abl). The simulation results showed that ossPTMetaD requires much fewer numbers of replicas than temperature REMD (T-REMD) with a reduction of ∼70% to achieve a similar exchange ratio. Although it does not obey the detailed balance condition, ossPTMetaD provides consistent results with T-REMD and experimental data. The high accessibility of the large conformational change of protein systems by ossPTMetaD, especially in simulating the very challenging DFG-motif flip of protein kinases, demonstrated its high efficiency and robustness in the characterization of the large-scale protein conformational change pathway and associated free energy profile.

All-atom adaptively biased path optimization of Src kinase conformational inactivation: Switched electrostatic network in the concerted motion of αC helix and the activation loop

A method to optimize a conformational pathway through a space of well-chosen reduced variables is employed to advance our understanding of protein conformational equilibrium. The adaptively biased path optimization strategy utilizes unrestricted, enhanced sampling in the region of a path in the reduced-variable space to identify a broad path between two stable end-states. Application to the inactivation transition of the Src tyrosine kinase catalytic domain reveals new insight into this well studied conformational equilibrium. The mechanistic description gained from identifying the motions and structural features along the path includes details of the switched electrostatic network found to underpin the transition. The free energy barrier along the path results from rotation of a helix, αC, that is tightly correlated with motions in the activation loop (A-loop) as well as distal regions in the C-lobe. Path profiles of the reduced variables clearly demonstrate the strongly correlated motions. The exchange of electrostatic interactions among residues in the network is key to these interdependent motions. In addition, the increased resolution from an all-atom model in defining the path shows multiple components for the A-loop motion and that different parts of the A-loop contribute throughout the length of the path.

Activation and Inactivation of the FLT3 Kinase: Pathway Intermediates and the Free Energy of Transition

The aberrant expression of kinases is often associated with pathologies such as cancer and autoimmune diseases. Like other types of enzymes, kinases can adopt active and inactive states, where a shift toward more stable active state often leads to disease. Dozens of kinase inhibitors are, therefore, used as drugs. Most of these bind to either the inactive or active state. In this work, we study the transitions between these two states in FLT3, an important drug target in leukemias. Kinases are composed of two lobes (N- and C-terminal lobes) with the catalytic site in-between. Through projection of the largest motions obtained through molecular dynamics (MD) simulations, we show that each of the end-states (active or inactive) already possess the ability for transition as the two lobes rotate which initiates the transition. A targeted simulation approach known as essential dynamics sampling (EDS) was used to speed up the transition between the two protein states. Coupling the EDS to implicit-solvent MD was performed to estimate the free energy barriers of the transitions. The activation energies were found in good agreement with previous estimates obtained for other kinases. Finally, we identified FLT3 intermediates that assumed configurations that resemble that of the c-Src nonreceptor tyrosine kinase. The intermediates show better binding to the drug ponatinib than c-Src and the inactive state of FLT3. This suggests that targeting intermediate states can be used to explain the drug-binding patterns of kinases and for rational drug design.

Direct cysteine sulfenylation drives activation of the Src kinase

The Src kinase controls aspects of cell biology and its activity is regulated by intramolecular structural changes induced by protein interactions and tyrosine phosphorylation. Recent studies indicate that Src is additionally regulated by redox-dependent mechanisms, involving oxidative modification(s) of cysteines within the Src protein, although the nature and molecular-level impact of Src cysteine oxidation are unknown. Using a combination of biochemical and cell-based studies, we establish the critical importance of two Src cysteine residues, Cys-185 and Cys-277, as targets for H2O2-mediated sulfenylation (Cys-SOH) in redox-dependent kinase activation in response to NADPH oxidase-dependent signaling. Molecular dynamics and metadynamics simulations reveal the structural impact of sulfenylation of these cysteines, indicating that Cys-277-SOH enables solvent exposure of Tyr-416 to promote its (auto)phosphorylation, and that Cys-185-SOH destabilizes pTyr-527 binding to the SH2 domain. These redox-dependent Src activation mechanisms offer opportunities for development of Src-selective inhibitors in treatment of diseases where Src is aberrantly activated.

Exploring Configuration Space and Path Space of Biomolecules Using Enhanced Sampling Techniques-Searching for Mechanism and Kinetics of Biomolecular Functions

To understand functions of biomolecules such as proteins, not only structures but their conformational change and kinetics need to be characterized, but its atomistic details are hard to obtain both experimentally and computationally. Here, we review our recent computational studies using novel enhanced sampling techniques for conformational sampling of biomolecules and calculations of their kinetics. For efficiently characterizing the free energy landscape of a biomolecule, we introduce the multiscale enhanced sampling method, which uses a combined system of atomistic and coarse-grained models. Based on the idea of Hamiltonian replica exchange, we can recover the statistical properties of the atomistic model without any biases. We next introduce the string method as a path search method to calculate the minimum free energy pathways along a multidimensional curve in high dimensional space. Finally we introduce novel methods to calculate kinetics of biomolecules based on the ideas of path sampling: one is the Onsager⁻Machlup action method, and the other is the weighted ensemble method. Some applications of the above methods to biomolecular systems are also discussed and illustrated.

The Src module: an ancient scaffold in the evolution of cytoplasmic tyrosine kinases

Tyrosine kinases were first discovered as the protein products of viral oncogenes. We now know that this large family of metazoan enzymes includes nearly one hundred structurally diverse members. Tyrosine kinases are broadly classified into two groups: the transmembrane receptor tyrosine kinases, which sense extracellular stimuli, and the cytoplasmic tyrosine kinases, which contain modular ligand-binding domains and propagate intracellular signals. Several families of cytoplasmic tyrosine kinases have in common a core architecture, the "Src module," composed of a Src-homology 3 (SH3) domain, a Src-homology 2 (SH2) domain, and a kinase domain. Each of these families is defined by additional elaborations on this core architecture. Structural, functional, and evolutionary studies have revealed a unifying set of principles underlying the activity and regulation of tyrosine kinases built on the Src module. The discovery of these conserved properties has shaped our knowledge of the workings of protein kinases in general, and it has had important implications for our understanding of kinase dysregulation in disease and the development of effective kinase-targeted therapies.

New insights into the structural dynamics of the kinase JNK3

In this work, we study the dynamics and the energetics of the all-atom structure of a neuronal-specific serine/threonine kinase c-Jun N-terminal kinase 3 (JNK3) in three states: unphosphorylated, phosphorylated, and ATP-bound phosphorylated. A series of 2 µs atomistic simulations followed by a conformational landscape mapping and a principal component analysis supports the mechanistic understanding of the JNK3 inactivation/activation process and also indicates key structural intermediates. Our analysis reveals that the unphosphorylated JNK3 undergoes the ‘open-to-closed’ movement via a two-step mechanism. Furthermore, the phosphorylation and ATP-binding allow the JNK3 kinase to attain a fully active conformation. JNK3 is a widely studied target for small-drugs used to treat a variety of neurological disorders. We believe that the mechanistic understanding of the large-conformational changes upon the activation of JNK3 will aid the development of novel targeted therapeutics.

Energetics and conformational pathways of functional rotation in the multidrug transporter AcrB

The multidrug transporter AcrB transports a broad range of drugs out of the cell by means of the proton-motive force. The asymmetric crystal structure of trimeric AcrB suggests a functionally rotating mechanism for drug transport. Despite various supportive forms of evidence from biochemical and simulation studies for this mechanism, the link between the functional rotation and proton translocation across the membrane remains elusive. Here, calculating the minimum free energy pathway of the functional rotation for the complete AcrB trimer, we describe the structural and energetic basis behind the coupling between the functional rotation and the proton translocation at atomic resolution. Free energy calculations show that protonation of Asp408 in the transmembrane portion of the drug-bound protomer drives the functional rotation. The conformational pathway identifies vertical shear motions among several transmembrane helices, which regulate alternate access of water in the transmembrane as well as peristaltic motions that pump drugs in the periplasm.

Millisecond dynamics of BTK reveal kinome-wide conformational plasticity within the apo kinase domain

Bruton tyrosine kinase (BTK) is a key enzyme in B-cell development whose improper regulation causes severe immunodeficiency diseases. Design of selective BTK therapeutics would benefit from improved, in-silico structural modeling of the kinase’s solution ensemble. However, this remains challenging due to the immense computational cost of sampling events on biological timescales. In this work, we combine multi-millisecond molecular dynamics (MD) simulations with Markov state models (MSMs) to report on the thermodynamics, kinetics, and accessible states of BTK’s kinase domain. Our conformational landscape links the active state to several inactive states, connected via a structurally diverse intermediate. Our calculations predict a kinome-wide conformational plasticity, and indicate the presence of several new potentially druggable BTK states. We further find that the population of these states and the kinetics of their inter-conversion are modulated by protonation of an aspartate residue, establishing the power of MD & MSMs in predicting effects of chemical perturbations.

Computational membrane biophysics: From ion channel interactions with drugs to cellular function

The rapid development of experimental and computational techniques has changed fundamentally our understanding of cellular-membrane transport. The advent of powerful computers and refined force-fields for proteins, ions, and lipids has expanded the applicability of Molecular Dynamics (MD) simulations. A myriad of cellular responses is modulated through the binding of endogenous and exogenous ligands (e.g. neurotransmitters and drugs, respectively) to ion channels. Deciphering the thermodynamics and kinetics of the ligand binding processes to these membrane proteins is at the heart of modern drug development. The ever-increasing computational power has already provided insightful data on the thermodynamics and kinetics of drug-target interactions, free energies of solvation, and partitioning into lipid bilayers for drugs. This review aims to provide a brief summary about modeling approaches to map out crucial binding pathways with intermediate conformations and free-energy surfaces for drug-ion channel binding mechanisms that are responsible for multiple effects on cellular functions. We will discuss post-processing analysis of simulation-generated data, which are then transformed to kinetic models to better understand the molecular underpinning of the experimental observables under the influence of drugs or mutations in ion channels. This review highlights crucial mathematical frameworks and perspectives on bridging different well-established computational techniques to connect the dynamics and timescales from all-atom MD and free energy simulations of ion channels to the physiology of action potentials in cellular models. This article is part of a Special Issue entitled: Biophysics in Canada, edited by Lewis Kay, John Baenziger, Albert Berghuis and Peter Tieleman.

Molecular determinants for the thermodynamic and functional divergence of uniporter GLUT1 and proton symporter XylE

GLUT1 facilitates the down-gradient translocation of D-glucose across cell membrane in mammals. XylE, an Escherichia coli homolog of GLUT1, utilizes proton gradient as an energy source to drive uphill D-xylose transport. Previous studies of XylE and GLUT1 suggest that the variation between an acidic residue (Asp27 in XylE) and a neutral one (Asn29 in GLUT1) is a key element for their mechanistic divergence. In this work, we combined computational and biochemical approaches to investigate the mechanism of proton coupling by XylE and the functional divergence between GLUT1 and XylE. Using molecular dynamics simulations, we evaluated the free energy profiles of the transition between inward- and outward-facing conformations for the apo proteins. Our results revealed the correlation between the protonation state and conformational preference in XylE, which is supported by the crystal structures. In addition, our simulations suggested a thermodynamic difference between XylE and GLUT1 that cannot be explained by the single residue variation at the protonation site. To understand the molecular basis, we applied Bayesian network models to analyze the alteration in the architecture of the hydrogen bond networks during conformational transition. The models and subsequent experimental validation suggest that multiple residue substitutions are required to produce the thermodynamic and functional distinction between XylE and GLUT1. Despite the lack of simulation studies with substrates, these computational and biochemical characterizations provide unprecedented insight into the mechanistic difference between proton symporters and uniporters.

We seek to address one intriguing question, the mechanistic distinction between active proton-coupled symporters and passive uniporters that are related in evolution. Proton-coupled symporters harness the transmembrane proton gradient to drive the substrate transport, while uniporters can only facilitate the passive substrate translocation. In this work, we focus on two sugar transporters GLUT1 and XylE, which belong to symporters and uniporters respectively but have high sequence similarity. We first applied molecular dynamics simulations to characterize the thermodynamic behaviors of apo GLUT1 and XylE, which are supposed to provide prominent details of mechanisms. From the identified difference in thermodynamics, we concluded that neutralizing protonation site in XylE is insufficient for its conversion to GLUT1 analog. To pinpoint extra elements contributing to their evolutionary divergence, we developed a novel network modeling scheme based on Bayesian network which shows impressive predictive power on residue mutations. Our models suggested the detailed mechanism of proton coupling in XylE and molecular basis of symporter/uniporter discrepancy. Furthermore, the modeling scheme could help to guide the design of biomolecules for desired functions.

Tyrosine Kinase Activation and Conformational Flexibility: Lessons from Src-Family Tyrosine Kinases

Protein kinases are enzymes that catalyze the covalent transfer of the γ-phosphate of an adenosine triphosphate (ATP) molecule onto a tyrosine, serine, threonine, or histidine residue in the substrate and thus send a chemical signal to networks of downstream proteins. They are important cellular signaling enzymes that regulate cell growth, proliferation, metabolism, differentiation, and migration. Unregulated protein kinase activity is often associated with a wide range of diseases, therefore making protein kinases major therapeutic targets. A prototypical system of central interest to understand the regulation of kinase activity is provided by tyrosine kinase c-Src, which belongs to the family of Src-related non-receptor tyrosine kinases (SFKs). Although the broad picture of autoinhibition via the regulatory domains and via the phosphorylation of the C-terminal tail is well characterized from a structural point of view, a detailed mechanistic understanding at the atomic-level is lacking. Advanced computational methods based on all-atom molecular dynamics (MD) simulations are employed to advance our understanding of tyrosine kinase activation. The computational studies suggest that the isolated kinase domain (KD) is energetically most favorable in the inactive conformation when the activation loop (A-loop) of the KD is not phosphorylated. The KD makes transient visits to a catalytically competent active-like conformation. The process of bimolecular trans-autophosphorylation of the A-loop eventually locks the KD in the active state. Activating point mutations may act by slightly increasing the population of the active-like conformation, enhancing the availability of the A-loop to be phosphorylated. The Src-homology 2 (SH2) and Src-homology 3 (SH3) regulatory domains, depending upon their configuration, either promote the inactive or the active state of the kinase domain. In addition to the roles played by the SH3, SH2, and KD, the Src-homology 4-Unique domain (SH4-U) region also serves as a key moderator of substrate specificity and kinase function. Thus, a fundamental understanding of the conformational propensity of the SH4-U region and how this affects the association to the membrane surface are likely to lead to the discovery of new intermediate states and alternate strategies for inhibition of kinase activity for drug discovery. The existence of a multitude of KD conformations poses a great challenge aimed at the design of specific inhibitors. One promising computational strategy to explore the conformational flexibility of the KD is to construct Markov state models from aggregated MD data.

The Activation of c-Src Tyrosine Kinase: Conformational Transition Pathway and Free Energy Landscape

Tyrosine kinases are important cellular signaling allosteric enzymes that regulate cell growth, proliferation, metabolism, differentiation, and migration. Their activity must be tightly controlled, and malfunction can lead to a variety of diseases, particularly cancer. The nonreceptor tyrosine kinase c-Src, a prototypical model system and a representative member of the Src-family, functions as complex multidomain allosteric molecular switches comprising SH2 and SH3 domains modulating the activity of the catalytic domain. The broad picture of self-inhibition of c-Src via the SH2 and SH3 regulatory domains is well characterized from a structural point of view, but a detailed molecular mechanism understanding is nonetheless still lacking. Here, we use advanced computational methods based on all-atom molecular dynamics simulations with explicit solvent to advance our understanding of kinase activation. To elucidate the mechanism of regulation and self-inhibition, we have computed the pathway and the free energy landscapes for the "inactive-to-active" conformational transition of c-Src for different configurations of the SH2 and SH3 domains. Using the isolated c-Src catalytic domain as a baseline for comparison, it is observed that the SH2 and SH3 domains, depending upon their bound orientation, promote either the inactive or active state of the catalytic domain. The regulatory structural information from the SH2-SH3 tandem is allosterically transmitted via the N-terminal linker of the catalytic domain. Analysis of the conformational transition pathways also illustrates the importance of the conserved tryptophan 260 in activating c-Src, and reveals a series of concerted events during the activation process.

Conformational Transitions and Alternating-Access Mechanism in the Sarcoplasmic Reticulum Calcium Pump

Ion pumps are integral membrane proteins responsible for transporting ions against concentration gradients across biological membranes. Sarco/endoplasmic reticulum Ca2+-ATPase (SERCA), a member of the P-type ATPases family, transports two calcium ions per hydrolyzed ATP molecule via an "alternating-access" mechanism. High-resolution crystallographic structures provide invaluable insight on the structural mechanism of the ion pumping process. However, to understand the molecular details of how ATP hydrolysis is coupled to calcium transport, it is necessary to gain knowledge about the conformational transition pathways connecting the crystallographically resolved conformations. Large-scale transitions in SERCA occur at time-scales beyond the current reach of unbiased molecular dynamics simulations. Here, we overcome this challenge by employing the string method, which represents a transition pathway as a chainofstates linking two conformational endpoints. Using a multiscale methodology, we have determined all-atom transition pathways for three main conformational transitions responsible for the alternating-access mechanism. The present pathways provide a clear chronology and ordering of the key events underlying the active transport of calcium ions by SERCA. Important conclusions are that the conformational transition that leads to occlusion with bound ATP and calcium is highly concerted and cooperative, the phosphorylation of Asp351 causes areorganization of the cytoplasmic domains that subsequently drives the opening of the luminal gate, and thereclosing of luminal gate induces a shift in the cytoplasmic domains that subsequently enables the dephosphorylation of Asp351-P. Formation of transient residue-residue contacts along the conformational transitions predicted by the computations provide an experimental route to test the general validity of the computational pathways.

Theoretical modeling of multiprotein complexes by iSPOT: Integration of small-angle X-ray scattering, hydroxyl radical footprinting, and computational docking

Structural determination of protein-protein complexes such as multidomain nuclear receptors has been challenging for high-resolution structural techniques. Here, we present a combined use of multiple biophysical methods, termed iSPOT, an integration of shape information from small-angle X-ray scattering (SAXS), protection factors probed by hydroxyl radical footprinting, and a large series of computationally docked conformations from rigid-body or molecular dynamics (MD) simulations. Specifically tested on two model systems, the power of iSPOT is demonstrated to accurately predict the structures of a large protein-protein complex (TGFβ-FKBP12) and a multidomain nuclear receptor homodimer (HNF-4α), based on the structures of individual components of the complexes. Although neither SAXS nor footprinting alone can yield an unambiguous picture for each complex, the combination of both, seamlessly integrated in iSPOT, narrows down the best-fit structures that are about 3.2Å and 4.2Å in RMSD from their corresponding crystal structures, respectively. Furthermore, this proof-of-principle study based on the data synthetically derived from available crystal structures shows that the iSPOT-using either rigid-body or MD-based flexible docking-is capable of overcoming the shortcomings of standalone computational methods, especially for HNF-4α. By taking advantage of the integration of SAXS-based shape information and footprinting-based protection/accessibility as well as computational docking, this iSPOT platform is set to be a powerful approach towards accurate integrated modeling of many challenging multiprotein complexes.

Microscopic Characterization of Membrane Transporter Function by In Silico Modeling and Simulation

Membrane transporters mediate one of the most fundamental processes in biology. They are the main gatekeepers controlling active traffic of materials in a highly selective and regulated manner between different cellular compartments demarcated by biological membranes. At the heart of the mechanism of membrane transporters lie protein conformational changes of diverse forms and magnitudes, which closely mediate critical aspects of the transport process, most importantly the coordinated motions of remotely located gating elements and their tight coupling to chemical processes such as binding, unbinding and translocation of transported substrate and cotransported ions, ATP binding and hydrolysis, and other molecular events fueling uphill transport of the cargo. An increasing number of functional studies have established the active participation of lipids and other components of biological membranes in the function of transporters and other membrane proteins, often acting as major signaling and regulating elements. Understanding the mechanistic details of these molecular processes require methods that offer high spatial and temporal resolutions. Computational modeling and simulations technologies empowered by advanced sampling and free energy calculations have reached a sufficiently mature state to become an indispensable component of mechanistic studies of membrane transporters in their natural environment of the membrane. In this article, we provide an overview of a number of major computational protocols and techniques commonly used in membrane transporter modeling and simulation studies. The article also includes practical hints on effective use of these methods, critical perspectives on their strengths and weak points, and examples of their successful applications to membrane transporters, selected from the research performed in our own laboratory.

Principles and Overview of Sampling Methods for Modeling Macromolecular Structure and Dynamics

Investigation of macromolecular structure and dynamics is fundamental to understanding how macromolecules carry out their functions in the cell. Significant advances have been made toward this end in silico, with a growing number of computational methods proposed yearly to study and simulate various aspects of macromolecular structure and dynamics. This review aims to provide an overview of recent advances, focusing primarily on methods proposed for exploring the structure space of macromolecules in isolation and in assemblies for the purpose of characterizing equilibrium structure and dynamics. In addition to surveying recent applications that showcase current capabilities of computational methods, this review highlights state-of-the-art algorithmic techniques proposed to overcome challenges posed in silico by the disparate spatial and time scales accessed by dynamic macromolecules. This review is not meant to be exhaustive, as such an endeavor is impossible, but rather aims to balance breadth and depth of strategies for modeling macromolecular structure and dynamics for a broad audience of novices and experts.

Evaluating the dynamics and electrostatic interactions of folded proteins in implicit solvents

Three implicit solvent models, namely GBMVII, FACTS, and SCPISM, were evaluated for their abilities to emulate an explicit solvent environment by comparing the simulated conformational ensembles, dynamics, and electrostatic interactions of the Src SH2 domain and the Lyn kinase domain. This assessment in terms of structural features in folded proteins expands upon the use of hydration energy as a metric for comparison. All-against-all rms coordinate deviation, average positional fluctuations, and ion-pair distance distribution were used to compare the implicit solvent models with the TIP3P explicit solvent model. Our study shows that the Src SH2 domains solvated with TIP3P, GBMVII, and FACTS sample similar global conformations. Additionally, the Src SH2 ion-pair distance distributions of solvent-exposed side chains corresponding to TIP3P, GBMVII, and FACTS do not differ substantially, indicating that GBMVII and FACTS are capable of modeling these electrostatic interactions. The ion-pair distance distributions of SCPISM are distinct from others, demonstrating that these electrostatic interactions are not adequately reproduced with the SCPISM model. On the other hand, for the Lyn kinase domain, a non-globular protein with bilobal structure and a large concavity on the surface, implicit solvent does not accurately model solvation to faithfully reproduce partially buried electrostatic interactions and lobe-lobe conformations. Our work reveals that local structure and dynamics of small, globular proteins are modeled well using FACTS and GBMVII. Nonetheless, global conformations and electrostatic interactions in concavities of multi-lobal proteins resulting from simulations with implicit solvent models do not match those obtained from explicit water simulations.

A stochastic roadmap method to model protein structural transitions

Evidence is emerging that the role of protein structure in disease needs to be rethought. Sequence mutations in proteins are often found to affect the rate at which a protein switches between structures. Modeling structural transitions in wildtype and variant proteins is central to understanding the molecular basis of disease. This paper investigates an efficient algorithmic realization of the stochastic roadmap simulation framework to model structural transitions in wildtype and variants of proteins implicated in human disorders. Our results indicate that the algorithm is able to extract useful information on the impact of mutations on protein structure and function.

Computational study of the W260A activating mutant of Src tyrosine kinase

Tyrosine kinases are enzymes playing a critical role in cellular signaling. Mutations causing increased in kinase activity are often associated with cancer and various pathologies. One example in Src tyrosine kinases is offered by the substitution of the highly conserved tryptophan 260 by an alanine (W260A), which has been shown to cause an increase in activity. Here, molecular dynamics simulations based on atomic models are carried out to characterize the conformational changes in the linker region and the catalytic (kinase) domain of Src kinase to elucidate the impact of the W260A mutation. Umbrella sampling calculations show that the conformation of the linker observed in the assembled down-regulated state of the kinase is most favored when the kinase domain is in the inactive state, whereas the conformation of the linker observed in the re-assembled up-regulated state of the kinase is favored when the kinase domain is in the unphosphorylated active-like state. The calculations further indicate that there are only small differences between the WT and W260A mutant. In both cases, the intermediates states are very similar and the down-regulated inactive conformation is the most stable state. However, the calculations also show that the free energy cost to reach the unphosphorylated active-like conformation is slightly smaller for the W260A mutant compared with WT. A simple kinetic model is developed and submitted to a Bayesian Monte Carlo analysis to illustrate how such small differences can contribute to accelerate the trans-autophosphorylation reaction and yield a large increase in the activity of the mutant as observed experimentally.

Theoretical and computational studies of peptides and receptors of the insulin family

Synergistic interactions among peptides and receptors of the insulin family are required for glucose homeostasis, normal cellular growth and development, proliferation, differentiation and other metabolic processes. The peptides of the insulin family are disulfide-linked single or dual-chain proteins, while receptors are ligand-activated transmembrane glycoproteins of the receptor tyrosine kinase (RTK) superfamily. Binding of ligands to the extracellular domains of receptors is known to initiate signaling via activation of intracellular kinase domains. While the structure of insulin has been known since 1969, recent decades have seen remarkable progress on the structural biology of apo and liganded receptor fragments. Here, we review how this useful structural information (on ligands and receptors) has enabled large-scale atomically-resolved simulations to elucidate the conformational dynamics of these biomolecules. Particularly, applications of molecular dynamics (MD) and Monte Carlo (MC) simulation methods are discussed in various contexts, including studies of isolated ligands, apo-receptors, ligand/receptor complexes and intracellular kinase domains. The review concludes with a brief overview and future outlook for modeling and computational studies in this family of proteins.

A dynamically coupled allosteric network underlies binding cooperativity in Src kinase

Protein tyrosine kinases are attractive drug targets because many human diseases are associated with the deregulation of kinase activity. However, how the catalytic kinase domain integrates different signals and switches from an active to an inactive conformation remains incompletely understood. Here we identify an allosteric network of dynamically coupled amino acids in Src kinase that connects regulatory sites to the ATP- and substrate-binding sites. Surprisingly, reactants (ATP and peptide substrates) bind with negative cooperativity to Src kinase while products (ADP and phosphopeptide) bind with positive cooperativity. We confirm the molecular details of the signal relay through the allosteric network by biochemical studies. Experiments on two additional protein tyrosine kinases indicate that the allosteric network may be largely conserved among these enzymes. Our work provides new insights into the regulation of protein tyrosine kinases and establishes a potential conduit by which resistance mutations to ATP-competitive kinase inhibitors can affect their activity.

Protein tyrosine kinases are subject to multiple regulatory mechanisms. Foda et al. show that reactants and products of the tyrosine kinase Src bind its catalytic domain with opposite cooperativity, and identify an allosteric network of dynamically coupled amino acids that underlie this behaviour.

Mapping the gating and permeation pathways in the voltage-gated proton channel Hv1

Voltage-gated proton channels (Hv1) are ubiquitous throughout nature and are implicated in numerous physiological processes. The gene encoding for Hv1, however, was only identified in 2006. The lack of sufficient structural information of this channel has hampered the understanding of the molecular mechanism of channel activation and proton permeation. This study uses both simulation and experimental approaches to further develop existing models of the Hv1 channel. Our study provides insights into features of channel gating and proton permeation pathway. We compare open- and closed-state structures developed previously with a recent crystal structure that traps the channel in a presumably closed state. Insights into gating pathways were provided using a combination of all-atom molecular dynamics simulations with a swarm of trajectories with the string method for extensive transition path sampling and evolution. A detailed residue-residue interaction profile and a hydration profile were studied to map the gating pathway in this channel. In particular, it allows us to identify potential intermediate states and compare them to the experimentally observed crystal structure of Takeshita et al. (Takeshita K, Sakata S, Yamashita E, Fujiwara Y, Kawanabe A, Kurokawa T, et al. X-ray crystal structure of voltage-gated proton channel. Nature 2014). The mechanisms governing ion transport in the wild-type and mutant Hv1 channels were studied by a combination of electrophysiological recordings and free energy simulations. With these results, we were able to further refine ideas about the location and function of the selectivity filter. The refined structural models will be essential for future investigations of this channel and the development of new drugs targeting cellular proton transport.

Enzymatic minimum free energy path calculations using swarms of trajectories

The development of approaches for simulating rare events in complex molecular systems is a central concern in chemical physics. In recent work, Roux and co-workers proposed a novel, swarms of trajectories (SoT) method for determining the transition paths of such events. It consists of the dynamical refinement on the system's free energy surface of a putative transition path that is parametrized in terms of a set of collective variables (CVs) that are identified as being important for the transition. In this work, we have implemented the SoT method and used it to investigate the catalytic mechanisms of two enzymatic reactions using hybrid QM/MM potentials. Our aim has been to test the performance of SoT for enzyme systems and to devise robust simulation protocols that can be employed in future studies of this type. We identify the conditions under which converged results can be obtained using inertial and Brownian dynamical evolutions of the CVs, show that the inclusion of several CVs can give significant additional insight into the mechanisms of the reactions, and show that the use of minimum energy paths as starting guesses can greatly accelerate path refinement.

Hotspot mutations in KIT receptor differentially modulate its allosterically coupled conformational dynamics: impact on activation and drug sensitivity

Receptor tyrosine kinase KIT controls many signal transduction pathways and represents a typical allosterically regulated protein. The mutation-induced deregulation of KIT activity impairs cellular physiological functions and causes serious human diseases. The impact of hotspots mutations (D816H/Y/N/V and V560G/D) localized in crucial regulatory segments, the juxtamembrane region (JMR) and the activation (A-) loop, on KIT internal dynamics was systematically studied by molecular dynamics simulations. The mutational outcomes predicted in silico were correlated with in vitro and in vivo activation rates and drug sensitivities of KIT mutants. The allosteric regulation of KIT in the native and mutated forms is described in terms of communication between the two remote segments, JMR and A-loop. A strong correlation between the communication profile and the structural and dynamical features of KIT in the native and mutated forms was established. Our results provide new insight on the determinants of receptor KIT constitutive activation by mutations and resistance of KIT mutants to inhibitors. Depiction of an intra-molecular component of the communication network constitutes a first step towards an integrated description of vast communication pathways established by KIT in physiopathological contexts.

Receptor tyrosine kinase KIT plays a crucial role in the regulation of cell signaling. This allosterically controlled activity may be affected by gain-of-function mutations that promote the development of several cancers. Identification of the molecular basis of KIT constitutive activation and allosteric regulation has inspired computational study of KIT hotspot mutations. In the present contribution, we investigated the mutation-induced effects on KIT conformational dynamics and intra-protein communication conditionally on the mutation location and the nature of the substituting amino acid. Our data elucidate that all studied mutations stabilize an inactive non-autoinhibited state of KIT over the inactive auto-inhibited state prevalent for the native protein. This shift in the protein conformational landscape promotes KIT constitutive activation. Our in silico analysis established correlations between the structural and dynamical effects induced by oncogenic mutations and the mutants auto-activation rates and drug sensitivities measured in vitro and in vivo. Particularly, the A-loop mutations stabilize the drug-resistant forms, while the JMR mutations may facilitate inhibitors binding to the active site. Cross-correlations established between local and long-range structural and dynamical effects demonstrate the allosteric character of the gain-of-function mutations mode of action.

Structural dynamic analysis of apo and ATP-bound IRAK4 kinase

Interleukin-1 receptor-associated kinases (IRAKs) are Ser/Thr protein kinases that play an important role as signaling mediators in the signal transduction facilitated by the Toll-like receptor (TLR) and interleukin-1 receptor families. Among IRAK family members, IRAK4 is one of the drug targets for diseases related to the TLR and IL-1R signaling pathways. Experimental evidence suggests that the IRAK4 kinase domain is phosphorylated in its activation loop at T342, T345, and S346 in the fully activated state. However, the molecular interactions of subdomains within the active and inactive IRAK4 kinase domain are poorly understood. Hence, we employed a long-range molecular dynamics (MD) simulation to compare apo IRAK4 kinase domains (phosphorylated and unphosphorylated) and ATP-bound phosphorylated IRAK4 kinase domains. The MD results strongly suggested that lobe uncoupling occurs in apo unphosphorylated IRAK4 kinase via the disruption of the R334/T345 and R310/T345 interaction. In addition, apo unphosphorylated trajectory result in high mobility, particularly in the N lobe, activation segment, helix αG, and its adjoining loops. The Asp-Phe-Gly (DFG) and His-Arg-Asp (HRD) conserved kinase motif analysis showed the importance of these motifs in IRAK4 kinase activation. This study provides important information on the structural dynamics of IRAK4 kinase, which will aid in inhibitor development.

Computational Recipe for Efficient Description of Large-Scale Conformational Changes in Biomolecular Systems

Characterizing large-scale structural transitions in biomolecular systems poses major technical challenges to both experimental and computational approaches. On the computational side, efficient sampling of the configuration space along the transition pathway remains the most daunting challenge. Recognizing this issue, we introduce a knowledge-based computational approach toward describing large-scale conformational transitions using (i) nonequilibrium, driven simulations combined with work measurements and (ii) free energy calculations using empirically optimized biasing protocols. The first part is based on designing mechanistically relevant, system-specific reaction coordinates whose usefulness and applicability in inducing the transition of interest are examined using knowledge-based, qualitative assessments along with nonequilirbrium work measurements which provide an empirical framework for optimizing the biasing protocol. The second part employs the optimized biasing protocol resulting from the first part to initiate free energy calculations and characterize the transition quantitatively. Using a biasing protocol fine-tuned to a particular transition not only improves the accuracy of the resulting free energies but also speeds up the convergence. The efficiency of the sampling will be assessed by employing dimensionality reduction techniques to help detect possible flaws and provide potential improvements in the design of the biasing protocol. Structural transition of a membrane transporter will be used as an example to illustrate the workings of the proposed approach.

Using multiscale preconditioning to accelerate the convergence of iterative molecular calculations

Iterative procedures for optimizing properties of molecular models often converge slowly owing to the computational cost of accurately representing features of interest. Here, we introduce a preconditioning scheme that allows one to use a less expensive model to guide exploration of the energy landscape of a more expensive model and thus speed the discovery of locally stable states of the latter. We illustrate our approach in the contexts of energy minimization and the string method for finding transition pathways. The relation of the method to other multilevel simulation techniques and possible extensions are discussed.

Exploring the conformational transitions of biomolecular systems using a simple two-state anisotropic network model

Biomolecular conformational transitions are essential to biological functions. Most experimental methods report on the long-lived functional states of biomolecules, but information about the transition pathways between these stable states is generally scarce. Such transitions involve short-lived conformational states that are difficult to detect experimentally. For this reason, computational methods are needed to produce plausible hypothetical transition pathways that can then be probed experimentally. Here we propose a simple and computationally efficient method, called ANMPathway, for constructing a physically reasonable pathway between two endpoints of a conformational transition. We adopt a coarse-grained representation of the protein and construct a two-state potential by combining two elastic network models (ENMs) representative of the experimental structures resolved for the endpoints. The two-state potential has a cusp hypersurface in the configuration space where the energies from both the ENMs are equal. We first search for the minimum energy structure on the cusp hypersurface and then treat it as the transition state. The continuous pathway is subsequently constructed by following the steepest descent energy minimization trajectories starting from the transition state on each side of the cusp hypersurface. Application to several systems of broad biological interest such as adenylate kinase, ATP-driven calcium pump SERCA, leucine transporter and glutamate transporter shows that ANMPathway yields results in good agreement with those from other similar methods and with data obtained from all-atom molecular dynamics simulations, in support of the utility of this simple and efficient approach. Notably the method provides experimentally testable predictions, including the formation of non-native contacts during the transition which we were able to detect in two of the systems we studied. An open-access web server has been created to deliver ANMPathway results.

Many biomolecules are like tiny molecular machines that need to change their shapes and visit many states to perform their biological functions. For a complete molecular understanding of a biological process, one needs to have information on the relevant stable states of the system in question, as well as the pathways by which the system travels from one state to another. We report here an efficient computational method that uses the knowledge of experimental structures of a pair of stable states in order to construct an energetically favoravle pathway between them. We adopt a simple representation of the molecular system by replacing the atoms with beads connected by springs and constructing an energy function with two minima around the end-states. We searched for the structure with highest energy that the system is most likely to visit during the transition and created two paths starting from this structure and proceeding toward the end-states. The combined result of these two paths is the minimum energy pathway between the two stable states. We apply this method to study important structural changes in one enzyme and three large proteins that transport small molecules and ions across the cell membrane.

Activation pathway of Src kinase reveals intermediate states as targets for drug design

Unregulated activation of Src kinases leads to aberrant signaling, uncontrolled growth, and differentiation of cancerous cells. Reaching a complete mechanistic understanding of large scale conformational transformations underlying the activation of kinases could greatly help in the development of therapeutic drugs for the treatment of these pathologies. In principle, the nature of conformational transition could be modeled in silico via atomistic molecular dynamics simulations, although this is very challenging due to the long activation timescales. Here, we employ a computational paradigm that couples transition pathway techniques and Markov state model-based massively distributed simulations for mapping the conformational landscape of c-src tyrosine kinase. The computations provide the thermodynamics and kinetics of kinase activation for the first time, and help identify key structural intermediates. Furthermore, the presence of a novel allosteric site in an intermediate state of c-src that could be potentially utilized for drug design is predicted.

Generalized Scalable Multiple Copy Algorithms for Molecular Dynamics Simulations in NAMD

Computational methodologies that couple the dynamical evolution of a set of replicated copies of a system of interest offer powerful and flexible approaches to characterize complex molecular processes. Such multiple copy algorithms (MCAs) can be used to enhance sampling, compute reversible work and free energies, as well as refine transition pathways. Widely used examples of MCAs include temperature and Hamiltonian-tempering replica-exchange molecular dynamics (T-REMD and H-REMD), alchemical free energy perturbation with lambda replica-exchange (FEP/λ-REMD), umbrella sampling with Hamiltonian replica exchange (US/H-REMD), and string method with swarms-of-trajectories conformational transition pathways. Here, we report a robust and general implementation of MCAs for molecular dynamics (MD) simulations in the highly scalable program NAMD built upon the parallel programming system Charm++. Multiple concurrent NAMD instances are launched with internal partitions of Charm++ and located continuously within a single communication world. Messages between NAMD instances are passed by low-level point-to-point communication functions, which are accessible through NAMD's Tcl scripting interface. The communication-enabled Tcl scripting provides a sustainable application interface for end users to realize generalized MCAs without modifying the source code. Illustrative applications of MCAs with fine-grained inter-copy communication structure, including global lambda exchange in FEP/λ-REMD, window swapping US/H-REMD in multidimensional order parameter space, and string method with swarms-of-trajectories were carried out on IBM Blue Gene/Q to demonstrate the versatility and massive scalability of the present implementation.

Comparison between Mean Forces and Swarms-of-Trajectories String Methods

The original formulation of the string method in collective variable space is compared with a recent variant called string method with swarms-of-trajectories. The assumptions made in the original method are revisited and the significance of the minimum free energy path (MFEP) is discussed in the context of reactive events. These assumptions are compared to those made in the string method with swarms-of-trajectories, and shown to be equivalent in a certain regime: in particular an expression for the path identified by the swarms-of-trajectories method is given and shown to be closely related to the MFEP. Finally, the algorithmic aspects of both methods are compared.