Weighted-ensemble Brownian dynamics simulations for protein association reactions

An ergodic theorem for the weighted ensemble method

We study weighted ensemble, an interacting particle method for sampling distributions of Markov chains that has been used in computational chemistry since the 1990s. Many important applications of weighted ensemble require the computation of long time averages. We establish the consistency of weighted ensemble in this setting by proving an ergodic theorem for time averages. As part of the proof, we derive explicit variance formulas that could be useful for optimizing the method.

Implications of the unfolded state in the folding energetics of heterogeneous-backbone protein mimetics

Sequence-encoded folding is the foundation of protein structure and is also possible in synthetic chains of artificial chemical composition. In natural proteins, the characteristics of the unfolded state are as important as those of the folded state in determining folding energetics. While much is known about folded structures adopted by artificial protein-like chains, corresponding information about the unfolded states of these molecules is lacking. Here, we report the consequences of altered backbone composition on the structure, stability, and dynamics of the folded and unfolded states of a compact helix-rich protein. Characterization through a combination of biophysical experiments and atomistic simulation reveals effects of backbone modification that depend on both the type of artificial monomers employed and where they are applied in sequence. In general, introducing artificial connectivity in a way that reinforces characteristics of the unfolded state ensemble of the prototype natural protein minimizes the impact of chemical changes on folded stability. These findings have implications in the design of protein mimetics and provide an atomically detailed picture of the unfolded state of a natural protein and artificial analogues under non-denaturing conditions.

Biophysical experiments and atomistic simulation reveal impacts of protein backbone alteration on the ensemble that defines the unfolded state. These effects have implications on folded stability of protein mimetics.

Creating Maps of the Ligand Binding Landscape for Kinetics-Based Drug Discovery

Simulations of ligand-protein interactions can be very useful for drug design and to gain biological insight. Full pathways of ligand-protein binding can be used to get information about ligand binding transition states, which form the rate-limiting step of the binding and release processes. However, these simulations are typically limited by the presence of large energy barriers that separate stable poses of interest. Here we describe a simulation protocol for exploring and analyzing landscapes of ligand-protein interactions that makes use of molecular docking, enhanced molecular simulation with the weighted ensemble algorithm, and network analysis. It can be accomplished using a modest cluster of graphics processing units and freely accessible software. This protocol focuses on the construction and analysis of a network model of ligand binding poses and provides links to resources that describe the other steps in more detail. The end result of this protocol is a map of the ligand-protein binding landscape that identifies transition states of the ligand binding pathway, as well as alternative bound poses that could be stabilized with modifications to the ligand.

Flow-regulated nucleation protrusion theory for collapsed polymers

The globular-stretch transition of a collapsed polymer in low strain rate elongational flow is studied using polymeric protrusion kinetics scaling laws and numerical simulation. Results demonstrate the influence of fluid flow on the occurrence probability of long-length thermally nucleated polymeric protrusions, which regulate collapsed polymer unfolding in low strain rate flows. Further, we estimate that the globular-stretch transition rate (k_{s}) in low strain rate (∈[over ̇]) elongational flows varies as k_{s}∼e^{-α∈[over ̇]^{-1}}. Results here reveal that the existing approach of neglecting the effects of fluid flow on thermally nucleated protrusions distribution is not valid for analyzing polymer unfolding behavior in low strain rate flows. Neglecting such an effect overestimates the constant α in the scaling law of transition rate (k_{s}∼e^{-α∈[over ̇]^{-1}}) by a factor of 2.

Markovian Weighted Ensemble Milestoning (M-WEM): Long-Time Kinetics from Short Trajectories

We introduce a rare-event sampling scheme, named Markovian Weighted Ensemble Milestoning (M-WEM), which inlays a weighted ensemble framework within a Markovian milestoning theory to efficiently calculate thermodynamic and kinetic properties of long-time-scale biomolecular processes from short atomistic molecular dynamics simulations. M-WEM is tested on the Müller-Brown potential model, the conformational switching in alanine dipeptide, and the millisecond time-scale protein-ligand unbinding in a trypsin-benzamidine complex. Not only can M-WEM predict the kinetics of these processes with quantitative accuracy but it also allows for a scheme to reconstruct a multidimensional free-energy landscape along additional degrees of freedom, which are not part of the milestoning progress coordinate. For the ligand-receptor system, the experimental residence time, association and dissociation kinetics, and binding free energy could be reproduced using M-WEM within a simulation time of a few hundreds of nanoseconds, which is a fraction of the computational cost of other currently available methods, and close to 4 orders of magnitude less than the experimental residence time. Due to the high accuracy and low computational cost, the M-WEM approach can find potential applications in kinetics and free-energy-based computational drug design.

Gaussian-Accelerated Molecular Dynamics with the Weighted Ensemble Method: A Hybrid Method Improves Thermodynamic and Kinetic Sampling

Gaussian-accelerated molecular dynamics (GaMD) is a well-established enhanced sampling method for molecular dynamics simulations that effectively samples the potential energy landscape of the system by adding a boost potential, which smoothens the surface and lowers the energy barriers between states. GaMD is unable to give time-dependent properties such as kinetics directly. On the other hand, the weighted ensemble (WE) method can efficiently sample transitions between states with its many weighted trajectories, which directly yield rates and pathways. However, convergence to equilibrium conditions remains a challenge for the WE method. Hence, we have developed a hybrid method that combines the two methods, wherein GaMD is first used to sample the potential energy landscape of the system and WE is subsequently used to further sample the potential energy landscape and kinetic properties of interest. We show that the hybrid method can sample both thermodynamic and kinetic properties more accurately and quickly compared to using either method alone.

Atomistic description of molecular binding processes based on returning probability theory

The efficiency of molecular binding such as host-guest binding is commonly evaluated in terms of kinetics, such as rate coefficients. In general, to compute the coefficient of the overall binding process, we need to consider both the diffusion of reactants and barrier crossing to reach the bound state. Here, we develop a methodology of quantifying the rate coefficient of binding based on molecular dynamics simulation and returning probability (RP) theory proposed by Kim and Lee [J. Chem. Phys. 131, 014503 (2009)]. RP theory provides a tractable formula of the rate coefficient in terms of the thermodynamic stability and kinetics of the intermediate state on a predefined reaction coordinate. In this study, the interaction energy between reactants is utilized as the reaction coordinate, enabling us to effectively describe the reactants' relative position and orientation on one-dimensional space. Application of this method to the host-guest binding systems, which consist of β-cyclodextrin and small guest molecules, yields the rate coefficients consistent with the experimental results.

#COVIDisAirborne: AI-Enabled Multiscale Computational Microscopy of Delta SARS-CoV-2 in a Respiratory Aerosol

We seek to completely revise current models of airborne transmission of respiratory viruses by providing never-before-seen atomic-level views of the SARS-CoV-2 virus within a respiratory aerosol. Our work dramatically extends the capabilities of multiscale computational microscopy to address the significant gaps that exist in current experimental methods, which are limited in their ability to interrogate aerosols at the atomic/molecular level and thus ob-scure our understanding of airborne transmission. We demonstrate how our integrated data-driven platform provides a new way of exploring the composition, structure, and dynamics of aerosols and aerosolized viruses, while driving simulation method development along several important axes. We present a series of initial scientific discoveries for the SARS-CoV-2 Delta variant, noting that the full scientific impact of this work has yet to be realized.

ACM REFERENCE FORMAT

Abigail Dommer 1† , Lorenzo Casalino 1† , Fiona Kearns 1† , Mia Rosenfeld 1 , Nicholas Wauer 1 , Surl-Hee Ahn 1 , John Russo, 2 Sofia Oliveira 3 , Clare Morris 1 , AnthonyBogetti 4 , AndaTrifan 5,6 , Alexander Brace 5,7 , TerraSztain 1,8 , Austin Clyde 5,7 , Heng Ma 5 , Chakra Chennubhotla 4 , Hyungro Lee 9 , Matteo Turilli 9 , Syma Khalid 10 , Teresa Tamayo-Mendoza 11 , Matthew Welborn 11 , Anders Christensen 11 , Daniel G. A. Smith 11 , Zhuoran Qiao 12 , Sai Krishna Sirumalla 11 , Michael O'Connor 11 , Frederick Manby 11 , Anima Anandkumar 12,13 , David Hardy 6 , James Phillips 6 , Abraham Stern 13 , Josh Romero 13 , David Clark 13 , Mitchell Dorrell 14 , Tom Maiden 14 , Lei Huang 15 , John McCalpin 15 , Christo- pherWoods 3 , Alan Gray 13 , MattWilliams 3 , Bryan Barker 16 , HarindaRajapaksha 16 , Richard Pitts 16 , Tom Gibbs 13 , John Stone 6 , Daniel Zuckerman 2 *, Adrian Mulholland 3 *, Thomas MillerIII 11,12 *, ShantenuJha 9 *, Arvind Ramanathan 5 *, Lillian Chong 4 *, Rommie Amaro 1 *. 2021. #COVIDisAirborne: AI-Enabled Multiscale Computational Microscopy ofDeltaSARS-CoV-2 in a Respiratory Aerosol. In Supercomputing '21: International Conference for High Perfor-mance Computing, Networking, Storage, and Analysis . ACM, New York, NY, USA, 14 pages. https://doi.org/finalDOI.

The prediction of protein-ligand unbinding for modern drug discovery

INTRODUCTION

Drug-target thermodynamic and kinetic information have perennially important roles in drug design. The prediction of protein-ligand unbinding, which can provide important kinetic information, in experiments continues to face great challenges. Uncovering protein-ligand unbinding through molecular dynamics simulations has become efficient and inexpensive with the progress and enhancement of computing power and sampling methods.

AREAS COVERED

In this review, various sampling methods for protein-ligand unbinding and their basic principles are firstly briefly introduced. Then, their applications in predicting aspects of protein-ligand unbinding, including unbinding pathways, dissociation rate constants, residence time and binding affinity, are discussed.

EXPERT OPINION

Although various sampling methods have been successfully applied in numerous systems, they still have shortcomings and deficiencies. Most enhanced sampling methods require researchers to possess a wealth of prior knowledge of collective variables or reaction coordinates. In addition, most systems studied at present are relatively simple, and the study of complex systems in real drug research remains greatly challenging. Through the combination of machine learning and enhanced sampling methods, prediction accuracy can be further improved, and some problems encountered in complex systems also may be solved.

A gentle introduction to the non-equilibrium physics of trajectories: Theory, algorithms, and biomolecular applications

Despite the importance of non-equilibrium statistical mechanics in modern physics and related fields, the topic is often omitted from undergraduate and core-graduate curricula. Key aspects of non-equilibrium physics, however, can be understood with a minimum of formalism based on a rigorous trajectory picture. The fundamental object is the ensemble of trajectories, a set of independent time-evolving systems, which easily can be visualized or simulated (e.g., for protein folding) and which can be analyzed rigorously in analogy to an ensemble of static system configurations. The trajectory picture provides a straightforward basis for understanding first-passage times, "mechanisms" in complex systems, and fundamental constraints on the apparent reversibility of complex processes. Trajectories make concrete the physics underlying the diffusion and Fokker-Planck partial differential equations. Last but not least, trajectory ensembles underpin some of the most important algorithms that have provided significant advances in biomolecular studies of protein conformational and binding processes.

Collective variable-based enhanced sampling and machine learning

ABSTRACT

Collective variable-based enhanced sampling methods have been widely used to study thermodynamic properties of complex systems. Efficiency and accuracy of these enhanced sampling methods are affected by two factors: constructing appropriate collective variables for enhanced sampling and generating accurate free energy surfaces. Recently, many machine learning techniques have been developed to improve the quality of collective variables and the accuracy of free energy surfaces. Although machine learning has achieved great successes in improving enhanced sampling methods, there are still many challenges and open questions. In this perspective, we shall review recent developments on integrating machine learning techniques and collective variable-based enhanced sampling approaches. We also discuss challenges and future research directions including generating kinetic information, exploring high-dimensional free energy surfaces, and efficiently sampling all-atom configurations.

GRAPHIC ABSTRACT

Toward Convergence in Free Energy Calculations for Protein Conformational Changes: A Case Study on the Thin Gate of Mhp1 Transporter

It has been challenging to obtain reliable free energies for protein conformational changes from all-atom molecular dynamics simulations, despite the availability of many enhanced sampling techniques. To alleviate the difficulties associated with the enormous complexity of the conformational space, here we propose a few practical strategies for such calculations, including (1) a stringent method to examine convergence by comparing independent simulations starting from different initial coordinates, (2) adoption of multistep schemes in which the complete conformational change consists of multiple transition steps, each sampled using a distinct reaction coordinate, and (3) application of boundary restraints to simplify the conformational space. We demonstrate these strategies on the conformational changes between the outward-facing and outward-occluded states of the Mhp1 membrane transporter, obtaining the equilibrium thermodynamics of the relevant metastable states, the kinetic rates between these states, and the reactive trajectories that reveal the atomic details of spontaneous transitions. Our approaches thus promise convergent and reliable calculations to examine intuition-based hypotheses and to eventually elucidate the underlying molecular mechanisms of reversible conformational changes in complex protein systems.

AI-driven prediction of SARS-CoV-2 variant binding trends from atomistic simulations

We present a novel technique to predict binding affinity trends between two molecules from atomistic molecular dynamics simulations. The technique uses a neural network algorithm applied to a series of images encoding the distance between two molecules in time. We demonstrate that our algorithm is capable of separating with high accuracy non-hydrophobic mutations with low binding affinity from those with high binding affinity. Moreover, we show high accuracy in prediction using a small subset of the simulation, therefore requiring a much shorter simulation time. We apply our algorithm to the binding between several variants of the SARS-CoV-2 spike protein and the human receptor ACE2.

AI-driven multiscale simulations illuminate mechanisms of SARS-CoV-2 spike dynamics

We develop a generalizable AI-driven workflow that leverages heterogeneous HPC resources to explore the time-dependent dynamics of molecular systems. We use this workflow to investigate the mechanisms of infectivity of the SARS-CoV-2 spike protein, the main viral infection machinery. Our workflow enables more efficient investigation of spike dynamics in a variety of complex environments, including within a complete SARS-CoV-2 viral envelope simulation, which contains 305 million atoms and shows strong scaling on ORNL Summit using NAMD. We present several novel scientific discoveries, including the elucidation of the spike's full glycan shield, the role of spike glycans in modulating the infectivity of the virus, and the characterization of the flexible interactions between the spike and the human ACE2 receptor. We also demonstrate how AI can accelerate conformational sampling across different systems and pave the way for the future application of such methods to additional studies in SARS-CoV-2 and other molecular systems.

A glycan gate controls opening of the SARS-CoV-2 spike protein

SARS-CoV-2 infection is controlled by the opening of the spike protein receptor binding domain (RBD), which transitions from a glycan-shielded 'down' to an exposed 'up' state to bind the human angiotensin-converting enzyme 2 receptor and infect cells. While snapshots of the 'up' and 'down' states have been obtained by cryo-electron microscopy and cryo-electron tomagraphy, details of the RBD-opening transition evade experimental characterization. Here over 130 µs of weighted ensemble simulations of the fully glycosylated spike ectodomain allow us to characterize more than 300 continuous, kinetically unbiased RBD-opening pathways. Together with ManifoldEM analysis of cryo-electron microscopy data and biolayer interferometry experiments, we reveal a gating role for the N-glycan at position N343, which facilitates RBD opening. Residues D405, R408 and D427 also participate. The atomic-level characterization of the glycosylated spike activation mechanism provided herein represents a landmark study for ensemble pathway simulations and offers a foundation for understanding the fundamental mechanisms of SARS-CoV-2 viral entry and infection.

The Next Frontier for Designing Switchable Proteins: Rational Enhancement of Kinetics

Designing proteins that can switch between active (ON) and inactive (OFF) conformations in response to signals such as ligand binding and incident light has been a tantalizing endeavor in protein engineering for over a decade. While such designs have yielded novel biosensors, therapeutic agents, and smart biomaterials, the response times (times for switching ON and OFF) of many switches have been too slow to be of practical use. Among the defining properties of such switches, the kinetics of switching has been the most challenging to optimize. This is largely due to the difficulty of characterizing the structures of transient states, which are required for manipulating the height of the effective free energy barrier between the ON and OFF states. We share our perspective of the most promising new experimental and computational strategies over the past several years for tackling this next frontier for designing switchable proteins.

Classification of protein-protein association rates based on biophysical informatics

BACKGROUND

Proteins form various complexes to carry out their versatile functions in cells. The dynamic properties of protein complex formation are mainly characterized by the association rates which measures how fast these complexes can be formed. It was experimentally observed that the association rates span an extremely wide range with over ten orders of magnitudes. Identification of association rates within this spectrum for specific protein complexes is therefore essential for us to understand their functional roles.

RESULTS

To tackle this problem, we integrate physics-based coarse-grained simulations into a neural-network-based classification model to estimate the range of association rates for protein complexes in a large-scale benchmark set. The cross-validation results show that, when an optimal threshold was selected, we can reach the best performance with specificity, precision, sensitivity and overall accuracy all higher than 70%. The quality of our cross-validation data has also been testified by further statistical analysis. Additionally, given an independent testing set, we can successfully predict the group of association rates for eight protein complexes out of ten. Finally, the analysis of failed cases suggests the future implementation of conformational dynamics into simulation can further improve model.

CONCLUSIONS

In summary, this study demonstrated that a new modeling framework that combines biophysical simulations with bioinformatics approaches is able to identify protein-protein interactions with low association rates from those with higher association rates. This method thereby can serve as a useful addition to a collection of existing experimental approaches that measure biomolecular recognition.

Reduced efficacy of a Src kinase inhibitor in crowded protein solution

The inside of a cell is highly crowded with proteins and other biomolecules. How proteins express their specific functions together with many off-target proteins in crowded cellular environments is largely unknown. Here, we investigate an inhibitor binding with c-Src kinase using atomistic molecular dynamics (MD) simulations in dilute as well as crowded protein solution. The populations of the inhibitor, 4-amino-5-(4-methylphenyl)-7-(t-butyl)pyrazolo[3,4-d]pyrimidine (PP1), in bulk solution and on the surface of c-Src kinase are reduced as the concentration of crowder bovine serum albumins (BSAs) increases. This observation is consistent with the reduced PP1 inhibitor efficacy in experimental c-Src kinase assays in addition with BSAs. The crowded environment changes the major binding pathway of PP1 toward c-Src kinase compared to that in dilute solution. This change is explained based on the population shift mechanism of local conformations near the inhibitor binding site in c-Src kinase.

Prediction of the Drug-Target Binding Kinetics for Flexible Proteins by Comparative Binding Energy Analysis

There is growing consensus that the optimization of the kinetic parameters for drug-protein binding leads to improved drug efficacy. Therefore, computational methods have been developed to predict kinetic rates and to derive quantitative structure-kinetic relationships (QSKRs). Many of these methods are based on crystal structures of ligand-protein complexes. However, a drawback is that each ligand-protein complex is usually treated as having a single structure. Here, we present a modification of COMparative BINding Energy (COMBINE) analysis, which uses the structures of ligand-protein complexes to predict binding parameters. We introduce the option of using multiple structures to describe each ligand-protein complex in COMBINE analysis and apply this to study the effects of protein flexibility on the derivation of dissociation rate constants (k_off) for inhibitors of p38 mitogen-activated protein (MAP) kinase, which has a flexible binding site. Multiple structures were obtained for each ligand-protein complex by performing docking to an ensemble of protein configurations obtained from molecular dynamics simulations. Coefficients to scale ligand-protein interaction energies determined from energy-minimized structures of ligand-protein complexes were obtained by partial least squares regression, and they allowed for the computation of k_off values. The QSKR model obtained using single, energy-minimized crystal structures for each ligand-protein complex had higher predictive power than the QSKR model obtained with multiple structures from ensemble docking. However, incorporation of ligand-protein flexibility helped to highlight additional ligand-protein interactions that lead to longer residence times, such as interactions with residues Arg67 and Asp168, which are close to the ligand in many crystal structures. These results show that COMBINE analysis is a promising method to guide the design of compounds that bind to flexible proteins with improved binding kinetics.

Free energy and kinetics of cAMP permeation through connexin26 via applied voltage and milestoning

The connexin family is a diverse group of highly regulated wide-pore channels permeable to biological signaling molecules. Despite the critical roles of connexins in mediating selective molecular signaling in health and disease, the basis of molecular permeation through these pores remains unclear. Here, we report the thermodynamics and kinetics of binding and transport of a second messenger, adenosine-3',5'-cyclophosphate (cAMP), through a connexin26 hemichannel (Cx26). First, inward and outward fluxes of cAMP molecules solvated in KCl solution were obtained from 4 μs of ± 200 mV simulations. These fluxes data yielded a single-channel permeability of cAMP and cAMP/K⁺ permeability ratio consistent with experimentally measured values. The results from voltage simulations were then compared with the potential of mean force (PMF) and the mean first passage times (MFPTs) of a single cAMP without voltage, obtained from a total of 16.5 μs of Voronoi-tessellated Markovian milestoning simulations. Both the voltage simulations and the milestoning simulations revealed two cAMP-binding sites, for which the binding constants K_D and dissociation rates k_off were computed from PMF and MFPTs. The protein dipole inside the pore produces an asymmetric PMF, reflected in unequal cAMP MFPTs in each direction once within the pore. The free energy profiles under opposite voltages were derived from the milestoning PMF and revealed the interplay between voltage and channel polarity on the total free energy. In addition, we show how these factors influence the cAMP dipole vector during permeation, and how cAMP affects the local and nonlocal pore diameter in a position-dependent manner.

Computational Modeling as a Tool to Investigate PPI: From Drug Design to Tissue Engineering

Protein-protein interactions (PPIs) mediate a large number of important regulatory pathways. Their modulation represents an important strategy for discovering novel therapeutic agents. However, the features of PPI binding surfaces make the use of structure-based drug discovery methods very challenging. Among the diverse approaches used in the literature to tackle the problem, linear peptides have demonstrated to be a suitable methodology to discover PPI disruptors. Unfortunately, the poor pharmacokinetic properties of linear peptides prevent their direct use as drugs. However, they can be used as models to design enzyme resistant analogs including, cyclic peptides, peptide surrogates or peptidomimetics. Small molecules have a narrower set of targets they can bind to, but the screening technology based on virtual docking is robust and well tested, adding to the computational tools used to disrupt PPI. We review computational approaches used to understand and modulate PPI and highlight applications in a few case studies involved in physiological processes such as cell growth, apoptosis and intercellular communication.

A glycan gate controls opening of the SARS-CoV-2 spike protein

SARS-CoV-2 infection is controlled by the opening of the spike protein receptor binding domain (RBD), which transitions from a glycan-shielded "down" to an exposed "up" state in order to bind the human ACE2 receptor and infect cells. While snapshots of the "up" and "down" states have been obtained by cryoEM and cryoET, details of the RBD opening transition evade experimental characterization. Here, over 130 μs of weighted ensemble (WE) simulations of the fully glycosylated spike ectodomain allow us to characterize more than 300 continuous, kinetically unbiased RBD opening pathways. Together with ManifoldEM analysis of cryo-EM data and biolayer interferometry experiments, we reveal a gating role for the N-glycan at position N343, which facilitates RBD opening. Residues D405, R408, and D427 also participate. The atomic-level characterization of the glycosylated spike activation mechanism provided herein achieves a new high-water mark for ensemble pathway simulations and offers a foundation for understanding the fundamental mechanisms of SARS-CoV-2 viral entry and infection.

What Markov State Models Can and Cannot Do: Correlation versus Path-Based Observables in Protein-Folding Models

Markov state models (MSMs) have been widely applied to study the kinetics and pathways of protein conformational dynamics based on statistical analysis of molecular dynamics (MD) simulations. These MSMs coarse-grain both configuration space and time in ways that limit what kinds of observables they can reproduce with high fidelity over different spatial and temporal resolutions. Despite their popularity, there is still limited understanding of which biophysical observables can be computed from these MSMs in a robust and unbiased manner, and which suffer from the space-time coarse-graining intrinsic in the MSM model. Most theoretical arguments and practical validity tests for MSMs rely on long-time equilibrium kinetics, such as the slowest relaxation time scales and experimentally observable time-correlation functions. Here, we perform an extensive assessment of the ability of well-validated protein folding MSMs to accurately reproduce path-based observable such as mean first-passage times (MFPTs) and transition path mechanisms compared to a direct trajectory analysis. We also assess a recently proposed class of history-augmented MSMs (haMSMs) that exploit additional information not accounted for in standard MSMs. We conclude with some practical guidance on the use of MSMs to study various problems in conformational dynamics of biomolecules. In brief, MSMs can accurately reproduce correlation functions slower than the lag time, but path-based observables can only be reliably reproduced if the lifetimes of states exceed the lag time, which is a much stricter requirement. Even in the presence of short-lived states, we find that haMSMs reproduce path-based observables more reliably.

Long-Time-Scale Predictions from Short-Trajectory Data: A Benchmark Analysis of the Trp-Cage Miniprotein

Elucidating physical mechanisms with statistical confidence from molecular dynamics simulations can be challenging owing to the many degrees of freedom that contribute to collective motions. To address this issue, we recently introduced a dynamical Galerkin approximation (DGA) [Thiede, E. H. J. Chem. Phys., 150, 2019, 244111], in which chemical kinetic statistics that satisfy equations of dynamical operators are represented by a basis expansion. Here, we reformulate this approach, clarifying (and reducing) the dependence on the choice of lag time. We present a new projection of the reactive current onto collective variables and provide improved estimators for rates and committors. We also present simple procedures for constructing suitable smoothly varying basis functions from arbitrary molecular features. To evaluate estimators and basis sets numerically, we generate and carefully validate a data set of short trajectories for the unfolding and folding of the trp-cage miniprotein, a well-studied system. Our analysis demonstrates a comprehensive strategy for characterizing reaction pathways quantitatively.

Path Ensembles for Pin1-Catalyzed Cis-Trans Isomerization of a Substrate Calculated by Weighted Ensemble Simulations

Pin1 enzyme protein recognizes specifically phosphorylated serine/threonine (pSer/pThr) and catalyzes the slow interconversion of the peptidyl-prolyl bond between cis and trans forms. Structural dynamics between the cis and trans forms are essential to reveal the underlying molecular mechanism of the catalysis. In this study, we apply the weighted ensemble (WE) simulation method to obtain comprehensive path ensembles for the Pin1-catalyzed isomerization process. Associated rate constants for both cis-to-trans and trans-to-cis isomerization are calculated to be submicroseconds time scales, which are in good agreement with the calculated free energy landscape where the cis form is slightly less favorable. The committor-like analysis indicates the shift of the transition state toward trans form (at the isomerization angle ω ∼ 110°) compared to the intrinsic position for the isolated substrate (ω ∼ 90°). The calculated structural ensemble clarifies a role of both the dual-histidine motif, His59/His157, and the basic residues, Lys63/Arg68/Arg69, to anchor both sides of the peptidyl-prolyl bond, the aromatic ring in Pro, and the phosphate in pSer, respectively. The rotation of the torsion angle is found to be facilitated by relaying the hydrogen-bond partner of the main-chain oxygen in pSer from Cys113 in the cis form to Arg68 in the trans form, through Ser154 at the transition state, which is really the cause of the shift in the transition state. The role of Ser154 as a driving force of the isomerization is confirmed by additional WE and free energy calculations for S154A mutant where the isomerization takes place slightly slower and the free energy barrier increases through the mutation. The present study shows the usefulness of the WE simulation for substantial path samplings between the reactant and product states, unraveling the molecular mechanism of the enzyme catalysis.

The RED scheme: Rate-constant estimation from pre-steady state weighted ensemble simulations

We present the Rate from Event Durations (RED) scheme, a new scheme that more efficiently calculates rate constants using the weighted ensemble path sampling strategy. This scheme enables rate-constant estimation from shorter trajectories by incorporating the probability distribution of event durations, or barrier-crossing times, from a simulation. We have applied the RED scheme to weighted ensemble simulations of a variety of rare-event processes that range in complexity: residue-level simulations of protein conformational switching, atomistic simulations of Na+/Cl- association in explicit solvent, and atomistic simulations of protein-protein association in explicit solvent. Rate constants were estimated with up to 50% greater efficiency than the original weighted ensemble scheme. Importantly, our scheme accounts for the systematic error that results from statistical bias toward the observation of events with short durations and reweights the event duration distribution accordingly. The RED scheme is relevant to any simulation strategy that involves unbiased trajectories of similar length to the most probable event duration, including weighted ensemble, milestoning, and standard simulations as well as the construction of Markov state models.

Predicting pathological von Willebrand factor unraveling in elongational flow

The globular-to-unraveled conformation transition of von Willebrand factor (vWF), a large polymeric glycoprotein in human blood plasma, is a crucial step in the process of clotting at sites of vascular injury. However, unraveling of vWF multimers in uninjured vasculature can lead to pathology (i.e., thrombus formation or degradation of vWF proteins by enzyme ADAMTS13, making them nonfunctional). To identify blood flow conditions that might induce pathological unraveling of vWF multimers, here we have computed the globular-to-unraveled transition rate of vWF multimers subjected to varying strain rate elongational flow by employing an enhanced sampling technique, the weighted ensemble method. Weighted ensemble sampling was employed instead of standard brute-force simulations because pathological blood flow conditions can induce undesired vWF unraveling on timescales potentially inaccessible to standard simulation methods. Results here indicate that brief but periodic exposure of vWF to the elongational flow of strain rate greater than or equal to 2500 s^-1 represents a source of possible pathology caused by the undesired unraveling of vWF multimers.

A Minimal, Adaptive Binning Scheme for Weighted Ensemble Simulations

A promising approach for simulating rare events with rigorous kinetics is the weighted ensemble path sampling strategy. One challenge of this strategy is the division of configurational space into bins for sampling. Here we present a minimal adaptive binning (MAB) scheme for the automated, adaptive placement of bins along a progress coordinate within the framework of the weighted ensemble strategy. Results reveal that the MAB binning scheme, despite its simplicity, is more efficient than a manual, fixed binning scheme in generating transitions over large free energy barriers, generating a diversity of pathways, estimating rate constants, and sampling conformations. The scheme is general and extensible to any rare-events sampling strategy that employs progress coordinates.

How μ-opioid receptor recognizes fentanyl

Roughly half of the drug overdose-related deaths in the United States are related to synthetic opioids represented by fentanyl which is a potent agonist of mu-opioid receptor (mOR). In recent years, X-ray crystal structures of mOR in complex with morphine derivatives have been determined; however, structural basis of mOR activation by fentanyl-like opioids remains lacking. Exploiting the X-ray structure of BU72-bound mOR and several molecular simulation techniques, we elucidated the detailed binding mechanism of fentanyl. Surprisingly, in addition to the salt-bridge binding mode common to morphinan opiates, fentanyl can move deeper and form a stable hydrogen bond with the conserved His2976.52, which has been suggested to modulate mOR's ligand affinity and pH dependence by previous mutagenesis experiments. Intriguingly, this secondary binding mode is only accessible when His2976.52 adopts a neutral HID tautomer. Alternative binding modes may represent a general mechanism in G protein-coupled receptor-ligand recognition.

An Implementation of Replica Exchange with Dynamical Scaling for Efficient Large-Scale Simulations

An implementation of the replica exchange with dynamical scaling (REDS) method in the commonly used molecular dynamics program GROMACS is presented. REDS is a replica exchange method that requires fewer replicas than conventional replica exchange while still providing data over a range of temperatures and can be used in either constant volume or constant pressure ensembles. Details for running REDS simulations are given, and an application to the human islet amyloid polypeptide (hIAPP) 11-25 fragment shows that the model efficiently samples conformational space.

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

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

Protein storytelling through physics

Every protein has a story-how it folds, what it binds, its biological actions, and how it misbehaves in aging or disease. Stories are often inferred from a protein's shape (i.e., its structure). But increasingly, stories are told using computational molecular physics (CMP). CMP is rooted in the principled physics of driving forces and reveals granular detail of conformational populations in space and time. Recent advances are accessing longer time scales, larger actions, and blind testing, enabling more of biology's stories to be told in the language of atomistic physics.

An Appraisal of Computer Simulation Approaches in Elucidating Biomolecular Recognition Pathways

Computer simulation approaches in biomolecular recognition processes have come a long way. In this Perspective, we highlight a series of recent success stories in which computer simulations have played a remarkable role in elucidating the atomic resolution mechanism of kinetic processes of protein-ligand binding in a quantitative fashion. In particular, we show that a robust combination of unbiased simulation, harnessed by a high-fidelity computing environment, and Markov state modeling approaches has been instrumental in revealing novel protein-ligand recognition pathways in multiple systems. We also elucidate the role of recent developments in enhanced sampling approaches in providing the much-needed impetus in accelerating simulation of the ligand recognition process. We identify multiple key issues, including force fields and the sampling bottleneck, which are currently preventing the field from achieving quantitative reconstruction of experimental measurements. Finally, we suggest a possible way forward via adoption of multiscale approaches and coarse-grained simulations as next steps toward efficient elucidation of ligand binding kinetics.

Molecular Dynamics Simulations in Drug Discovery and Pharmaceutical Development

Molecular dynamics (MD) simulations have become increasingly useful in the modern drug development process. In this review, we give a broad overview of the current application possibilities of MD in drug discovery and pharmaceutical development. Starting from the target validation step of the drug development process, we give several examples of how MD studies can give important insights into the dynamics and function of identified drug targets such as sirtuins, RAS proteins, or intrinsically disordered proteins. The role of MD in antibody design is also reviewed. In the lead discovery and lead optimization phases, MD facilitates the evaluation of the binding energetics and kinetics of the ligand-receptor interactions, therefore guiding the choice of the best candidate molecules for further development. The importance of considering the biological lipid bilayer environment in the MD simulations of membrane proteins is also discussed, using G-protein coupled receptors and ion channels as well as the drug-metabolizing cytochrome P450 enzymes as relevant examples. Lastly, we discuss the emerging role of MD simulations in facilitating the pharmaceutical formulation development of drugs and candidate drugs. Specifically, we look at how MD can be used in studying the crystalline and amorphous solids, the stability of amorphous drug or drug-polymer formulations, and drug solubility. Moreover, since nanoparticle drug formulations are of great interest in the field of drug delivery research, different applications of nano-particle simulations are also briefly summarized using multiple recent studies as examples. In the future, the role of MD simulations in facilitating the drug development process is likely to grow substantially with the increasing computer power and advancements in the development of force fields and enhanced MD methodologies.

Milestoning: An Efficient Approach for Atomically Detailed Simulations of Kinetics in Biophysics

Recent advances in theory and algorithms for atomically detailed simulations open the way to the study of the kinetics of a wide range of molecular processes in biophysics. The theories propose a shift from the traditionally very long molecular dynamic trajectories, which are exact but may not be efficient in the study of kinetics, to the use of a large number of short trajectories. The short trajectories exploit a mapping to a mesh in coarse space and allow for efficient calculations of kinetics and thermodynamics. In this review, I focus on one theory: Milestoning is a theory and an algorithm that offers a hierarchical calculation of properties of interest, such as the free energy profile and the mean first passage time. Approximations to the true long-time dynamics can be computed efficiently and assessed at different steps of the investigation. The theory is discussed and illustrated using two biophysical examples: ion permeation through a phospholipid membrane and protein translocation through a channel.

Wepy: A Flexible Software Framework for Simulating Rare Events with Weighted Ensemble Resampling

Here, we introduce the open-source software framework wepy (https://github.com/ADicksonLab/wepy) which is a toolkit for running and analyzing weighted ensemble (WE) simulations. The wepy toolkit is in pure Python and as such is highly portable and extensible, making it an excellent platform to develop and use new WE resampling algorithms such as WExplore, REVO, and others while leveraging the entire Python ecosystem. In addition, wepy simplifies WE-specific analyses by defining out-of-core tree-like data structures using the cross-platform HDF5 file format. In this paper, we discuss the motivations and challenges for simulating rare events in biomolecular systems. As has previously been shown, high-dimensional WE resampling algorithms such as WExplore and REVO have been successful at these tasks, especially for rare events that are difficult to describe by one or two collective variables. We explain in detail how wepy facilitates implementation of these algorithms, as well as aids in analyzing the unique structure of WE simulation results. To explain how wepy and WE work in general, we describe the mathematical formalism of WE, an overview of the architecture of wepy, and provide code examples of how to construct, run, and analyze simulation results for a protein-ligand system (T4 Lysozyme in an implicit solvent). This paper is written with a variety of readers in mind, including (1) those curious about how to leverage WE rare-event simulations for their domain, (2) current WE users who want to begin using new high-dimensional resamplers such as WExplore and REVO, and (3) expert users who would like to prototype or implement their own algorithms that can be easily adopted by others.

Thermodynamics and Kinetics of Drug-Target Binding by Molecular Simulation

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

How μ-Opioid Receptor Recognizes Fentanyl

In 2019, drug overdose has claimed over 70,000 lives in the United States. More than half of the deaths are related to synthetic opioids represented by fentanyl which is a potent agonist of mu-opioid receptor (mOR). In recent years, the crystal structures of mOR in complex with morphine derivatives have been determined; however, structural basis of mOR activation by fentanyl-like synthetic opioids remains lacking. Exploiting the X-ray structure of mOR bound to a morphinan ligand and several state-of-the-art simulation techniques, including weighted ensemble and continuous constant pH molecular dynamics, we elucidated the detailed binding mechanism of fentanyl with mOR. Surprisingly, in addition to forming a salt-bridge with Asp1473.32 in the orthosteric site common to morphinan opiates, fentanyl can move deeper and bind mOR through hydrogen bonding with a conserved histidine His2976.52, which has been shown to modulate mOR's ligand affinity and pH dependence in mutagenesis experiments, but its precise role remains unclear. Intriguingly, the secondary binding mode is only accessible when His297 adopts a neutral HID tautomer. Alternative binding modes and involvement of tautomer states may represent general mechanisms in G protein-coupled receptor (GPCR)-ligand recognition. Our work provides a starting point for understanding the molecular basis of mOR activation by fentanyl which has many analogs emerging at a rapid pace. The knowledge may also inform the design of safer analgesics to combat the opioid crisis. Current protein simulation studies employ standard protonation and tautomer states; our work demonstrates the need to move beyond the practice to advance our understanding of protein-ligand recognition.

Ranking of Ligand Binding Kinetics Using a Weighted Ensemble Approach and Comparison with a Multiscale Milestoning Approach

To improve lead optimization efforts in finding the right ligand, pharmaceutical industries need to know the ligand's binding kinetics, such as binding and unbinding rate constants, which often correlate with the ligand's efficacy in vivo. To predict binding kinetics efficiently, enhanced sampling methods, such as milestoning and the weighted ensemble (WE) method, have been used in molecular dynamics (MD) simulations of these systems. However, a comparison of these enhanced sampling methods in ranking ligands has not been done. Hence, a WE approach called the concurrent adaptive sampling (CAS) algorithm that uses MD simulations was used to rank seven ligands for β-cyclodextrin, a system in which a multiscale milestoning approach called simulation enabled estimation of kinetic rates (SEEKR) was also used, which uses both MD and Brownian dynamics simulations. Overall, the CAS algorithm can successfully rank ligands using the unbinding rate constant k_off values and binding free energy ΔG values, as SEEKR did, with reduced computational cost that is about the same as SEEKR. We compare the CAS algorithm simulations with different parameters and discuss the impact of parameters in ranking ligands and obtaining rate constant and binding free energy estimates. We also discuss similarities and differences and advantages and disadvantages of SEEKR and the CAS algorithm for future use.

AI-Driven Multiscale Simulations Illuminate Mechanisms of SARS-CoV-2 Spike Dynamics

We develop a generalizable AI-driven workflow that leverages heterogeneous HPC resources to explore the time-dependent dynamics of molecular systems. We use this workflow to investigate the mechanisms of infectivity of the SARS-CoV-2 spike protein, the main viral infection machinery. Our workflow enables more efficient investigation of spike dynamics in a variety of complex environments, including within a complete SARS-CoV-2 viral envelope simulation, which contains 305 million atoms and shows strong scaling on ORNL Summit using NAMD. We present several novel scientific discoveries, including the elucidation of the spike's full glycan shield, the role of spike glycans in modulating the infectivity of the virus, and the characterization of the flexible interactions between the spike and the human ACE2 receptor. We also demonstrate how AI can accelerate conformational sampling across different systems and pave the way for the future application of such methods to additional studies in SARS-CoV-2 and other molecular systems.

Accelerated Estimation of Long-Timescale Kinetics from Weighted Ensemble Simulation via Non-Markovian "Microbin" Analysis

The weighted ensemble (WE) simulation strategy provides unbiased sampling of nonequilibrium processes, such as molecular folding or binding, but the extraction of rate constants relies on characterizing steady-state behavior. Unfortunately, WE simulations of sufficiently complex systems will not relax to steady state on observed simulation times. Here, we show that a postsimulation clustering of molecular configurations into "microbins" using methods developed in the Markov State Model (MSM) community can yield unbiased kinetics from WE data before steady-state convergence of the WE simulation itself. Because WE trajectories are directional and not equilibrium distributed, the history-augmented MSM (haMSM) formulation can be used, which yields the mean first-passage time (MFPT) without bias for arbitrarily small lag times. Accurate kinetics can be obtained while bypassing the often prohibitive convergence requirements of the nonequilibrium weighted ensemble. We validate the method in a simple diffusive process on a two-dimensional (2D) random energy landscape and then analyze atomistic protein folding simulations using WE molecular dynamics. We report significant progress toward the unbiased estimation of protein folding times and pathways, though key challenges remain.

Kinetics and free energy of ligand dissociation using weighted ensemble milestoning

We consider the recently developed weighted ensemble milestoning (WEM) scheme [D. Ray and I. Andricioaei, J. Chem. Phys. 152, 234114 (2020)] and test its capability of simulating ligand-receptor dissociation dynamics. We performed WEM simulations on the following host-guest systems: Na⁺/Cl^- ion pair and 4-hydroxy-2-butanone ligand with FK506 binding protein. As a proof of principle, we show that the WEM formalism reproduces the Na⁺/Cl^- ion pair dissociation timescale and the free energy profile obtained from long conventional MD simulation. To increase the accuracy of WEM calculations applied to kinetics and thermodynamics in protein-ligand binding, we introduced a modified WEM scheme called weighted ensemble milestoning with restraint release (WEM-RR), which can increase the number of starting points per milestone without adding additional computational cost. WEM-RR calculations obtained a ligand residence time and binding free energy in agreement with experimental and previous computational results. Moreover, using the milestoning framework, the binding time and rate constants, dissociation constants, and committor probabilities could also be calculated at a low computational cost. We also present an analytical approach for estimating the association rate constant (k_on) when binding is primarily diffusion driven. We show that the WEM method can efficiently calculate multiple experimental observables describing ligand-receptor binding/unbinding and is a promising candidate for computer-aided inhibitor design.

Enhanced Jarzynski free energy calculations using weighted ensemble

The free energy of transitions between stable states is the key thermodynamic quantity that governs the relative probabilities of the forward and reverse reactions and the ratio of state probabilities at equilibrium. The binding free energy of a drug and its receptor is of particular interest, as it serves as an optimization function for drug design. Over the years, many computational methods have been developed to calculate binding free energies, and while many of these methods have a long history, issues such as convergence of free energy estimates and the projection of a binding process onto order parameters remain. Over 20 years ago, the Jarzynski equality was derived with the promise to calculate equilibrium free energies by measuring the work applied to short nonequilibrium trajectories. However, these calculations were found to be dominated by trajectories with low applied work that occur with extremely low probability. Here, we examine the combination of weighted ensemble algorithms with the Jarzynski equality. In this combined method, an ensemble of nonequilibrium trajectories are run in parallel, and cloning and merging operations are used to preferentially sample low-work trajectories that dominate the free energy calculations. Two additional methods are also examined: (i) a novel weighted ensemble resampler that samples trajectories directly according to their importance to the work of work and (ii) the diffusion Monte Carlo method using the applied work as the selection potential. We thoroughly examine both the accuracy and efficiency of unbinding free energy calculations for a series of model Lennard-Jones atom pairs with interaction strengths ranging from 2 kcal/mol to 20 kcal/mol. We find that weighted ensemble calculations can more efficiently determine accurate binding free energies, especially for deeper Lennard-Jones well depths.

Rare events and first passage time statistics from the energy landscape

We analyze the probability distribution of rare first passage times corresponding to transitions between product and reactant states in a kinetic transition network. The mean first passage times and the corresponding rate constants are analyzed in detail for two model landscapes and the double funnel landscape corresponding to an atomic cluster. Evaluation schemes based on eigendecomposition and kinetic path sampling, which both allow access to the first passage time distribution, are benchmarked against mean first passage times calculated using graph transformation. Numerical precision issues severely limit the useful temperature range for eigendecomposition, but kinetic path sampling is capable of extending the first passage time analysis to lower temperatures, where the kinetics of interest constitute rare events. We then investigate the influence of free energy based state regrouping schemes for the underlying network. Alternative formulations of the effective transition rates for a given regrouping are compared in detail to determine their numerical stability and capability to reproduce the true kinetics, including recent coarse-graining approaches that preserve occupancy cross correlation functions. We find that appropriate regrouping of states under the simplest local equilibrium approximation can provide reduced transition networks with useful accuracy at somewhat lower temperatures. Finally, a method is provided to systematically interpolate between the local equilibrium approximation and exact intergroup dynamics. Spectral analysis is applied to each grouping of states, employing a moment-based mode selection criterion to produce a reduced state space, which does not require any spectral gap to exist, but reduces to gap-based coarse graining as a special case. Implementations of the developed methods are freely available online.

Enhancing the promiscuity of a member of the Caspase protease family by rational design

The N-terminal cleavage of fusion tags to restore the native N-terminus of recombinant proteins is a challenging task and up to today, protocols need to be optimized for different proteins individually. Within this work, we present a novel protease that was designed in-silico to yield enhanced promiscuity toward different N-terminal amino acids. Two mutations in the active-site amino acids of human Caspase-2 were determined to increase the recognition of branched amino-acids, which show only poor binding capabilities in the unmutated protease. These mutations were determined by sequential and structural comparisons of Caspase-2 and Caspase-3 and their effect was additionally predicted using free-energy calculations. The two mutants proposed in the in-silico studies were expressed and in-vitro experiments confirmed the simulation results. Both mutants showed not only enhanced activities toward branched amino acids, but also smaller, unbranched amino acids. We believe that the created mutants constitute an important step toward generalized procedures to restore original N-termini of recombinant fusion proteins.

How mu-Opioid Receptor Recognizes Fentanyl

The opioid crisis has escalated during the COVID-19 pandemic. More than half of the overdose-related deaths are related to synthetic opioids represented by fentanyl which is a potent agonist of mu-opioid receptor (mOR). In recent years, crystal structures of mOR complexed with morphine derivatives have been determined; however, structural basis of mOR activation by fentanyl-like synthetic opioids remains lacking. Exploiting the X-ray structure of mOR bound to a morphinan ligand and several state-of-the-art simulation techniques, including weighted ensemble and continuous constant pH molecular dynamics, we elucidated the detailed binding mechanism of fentanyl with mOR. Surprisingly, in addition to the orthosteric site common to morphinan opiates, fentanyl can move deeper and bind mOR through hydrogen bonding with a conserved histidine H297, which has been shown to modulate mOR's ligand affinity and pH dependence in mutagenesis experiments, but its precise role remains unclear. Intriguingly, the secondary binding mode is only accessible when H297 adopts a neutral HID tautomer. Alternative binding modes and involvement of tautomer states may represent general mechanisms in G protein-coupled receptor (GPCR)-ligand recognition. Our work provides a starting point for understanding mOR activation by fentanyl analogs that are emerging at a rapid pace and assisting the design of safer analgesics to combat the opioid crisis. Current protein simulation studies employ standard protonation and tautomer states; our work demonstrates the need to move beyond the practice to advance our understanding of protein-ligand recognition.

Understanding the Impacts of Conformational Dynamics on the Regulation of Protein-Protein Association by a Multiscale Simulation Method

Complexes formed among diverse proteins carry out versatile functions in nearly all physiological processes. Association rates which measure how fast proteins form various complexes are of fundamental importance to characterize their functions. The association rates are not only determined by the energetic features at binding interfaces of a protein complex but also influenced by the intrinsic conformational dynamics of each protein in the complex. Unfortunately, how this conformational effect regulates protein association has never been calibrated on a systematic level. To tackle this problem, we developed a multiscale strategy to incorporate the information on protein conformational variations from Langevin dynamic simulations into a kinetic Monte Carlo algorithm of protein-protein association. By systematically testing this approach against a large-scale benchmark set, we found the association of a protein complex with a relatively rigid structure tends to be reduced by its conformational fluctuations. With specific examples, we further show that higher degrees of structural flexibility in various protein complexes can facilitate the searching and formation of intermolecular interactions and thereby accelerate their associations. In general, the integration of conformational dynamics can improve the correlation between experimentally measured association rates and computationally derived association probabilities. Finally, we analyzed the statistical distributions of different secondary structural types on protein-protein binding interfaces and their preference to the change of association rates. Our study, to the best of our knowledge, is the first computational method that systematically estimates the impacts of protein conformational dynamics on protein-protein association. It throws lights on the molecular mechanisms of how protein-protein recognition is kinetically modulated.