Water regulates the residence time of Benzamidine in Trypsin

Lambda-ABF-OPES: Faster Convergence with High Accuracy in Alchemical Free Energy Calculations

Predicting the binding affinity between small molecules and target macromolecules while combining both speed and accuracy is a cornerstone of modern computational drug discovery, which is critical for accelerating therapeutic development. Despite recent progress in molecular dynamics (MD) simulations, such as advanced polarizable force fields and enhanced sampling techniques, estimating absolute binding free energies (ABFEs) remains computationally challenging. To overcome these difficulties, we introduce a highly efficient hybrid methodology that couples the Lambda-adaptive biasing force (Lambda-ABF) scheme with on-the-fly probability enhanced sampling (OPES). This approach achieves up to a 9-fold improvement in sampling efficiency and computational speed compared to the original Lambda-ABF when used in conjunction with the AMOEBA polarizable force field, yielding converged results at a fraction of the cost of standard techniques.

Trust the force, but pull wisely: Structural insights into non-equilibrium response forces from pulling MD simulations

We analyze the coarse-grained equations of motion of molecular systems subject to external driving. As exemplary processes, we study by means of targeted and steered molecular dynamics simulations the dissociation of a sodium-chloride ion pair in water and ligand-protein unbinding of trypsin-benzamidine. We derive an exact generalization of Mori's Langevin equation that contains the memory kernel of the stationary process, an additive driving force, and a non-equilibrium response force describing the effects of the perturbed environment. We show that both the fluctuating force in the stationary case and the non-equilibrium response force in the driven cases exhibit spatial structure in their first and second moments. The latter depends sensitively on the employed driving protocols. For sodium chloride, we find that the first moment of the non-equilibrium response force matches the mean force for slow constrained pulling. In contrast, for all tested restrained pulling protocols, significant differences arise between the two properties in both systems. We conclude that the non-equilibrium response of the solvent needs to be taken into account carefully when analyzing data from pulling simulations.

Nonequilibrium Binding Free Energy Simulations: Minimizing Dissipation

While nowadays approaches for equilibrium free energy estimation are well established, nonequilibrium simulations represent both an appealing computational opportunity and a challenge. This kind of simulations allows for a trivially parallel scheme, but at the same time the significant amount of irreversible work often generated during the steering process (either alchemical or physical) can hinder the convergence of free energy estimators. Here, we discuss in depth this issue for the protein-ligand binding free energy estimation carried out via physical paths. We found that the water model and the parametrization of the path collective variables have a remarkable impact on the convergence rate of the estimators (e.g., Crooks). Finally, we provide practical recipes to enhance the convergence speed and minimize dissipation.

Long-Range Electrostatics in Serine Proteases: Machine Learning-Driven Reaction Sampling Yields Insights for Enzyme Design

Computational enzyme design is a promising technique for producing novel enzymes for industrial and clinical needs. A key challenge that this technique faces is to consistently achieve the desired activity. Fundamental studies of natural enzymes revealed critical contributions from second-shell - and even more distant - residues to their remarkable efficiency. In particular, such residues organize the internal electrostatic field to promote the reaction. Engineering such fields computationally proved to be a promising strategy, which, however, has some limitations. Charged residues necessarily form specific patterns of local interactions that may be exploited for structural integrity. As a result, it is impossible to probe the electrostatic field alone by substituting amino acids. We hypothesize that an approach that isolates the influences of residues' charges from other influences could yield deeper insights. We use molecular modeling with AI-enhanced QM/MM reaction sampling to implement such an approach and apply it to a model serine protease subtilisin. We find that the negative charge 8 Å away from the catalytic site is crucial to achieving the enzyme's catalytic efficiency, contributing more than 2 kcal/mol to lowering the barrier. In contrast, a positive charge from the second-closest charged residue opposes the efficiency of the reaction by raising the barrier by 0.8 kcal/mol. This result invites discussion into the role of this residue and trade-offs that might have taken place in the evolution of such enzymes. Our approach is transferable and can help investigate the evolution of electrostatic preorganization in other enzymes. We believe that the study and engineering of electrostatic fields in enzymes is a promising direction to advance both fundamental and applied enzymology and lead to the design of new powerful biocatalysts.

Enhanced Sampling Simulations of RNA-Peptide Binding Using Deep Learning Collective Variables

Enhanced sampling (ES) simulations of biomolecular recognition, such as binding small molecules to proteins and nucleic acid targets, protein-protein association, and protein-nucleic acid interactions, have gained significant attention in the simulation community because of their ability to sample long-time scale processes. However, a key challenge in implementing collective variable (CV)-based enhanced sampling methods is the selection of appropriate CVs that can distinguish the system's metastable states and, when biased, can effectively sample these states. This challenge is particularly acute when the binding of a flexible molecule to a conformationally rich host molecule is simulated, such as the binding of a peptide to an RNA. In such cases, a large number of CVs are required to capture the conformations of both the host and the guest as well as the binding process. Using such a large number of descriptors is impractical in any enhanced sampling simulation method. In our work, we used the recently developed deep targeted discriminant analysis (Deep-TDA) method to design CVs to study the binding of a cyclic peptide, L22, to a TAR RNA of HIV, which is a prototypical system. The Deep-TDA CV, obtained from a nonlinear combination of important contact pairs between the L22 peptide and the host RNA backbone atoms, along with the RNA apical loop RMSD as the second CV were used in the on-the-fly probability-based enhanced sampling (OPES) simulation to sample the reversible binding and unbinding of the L22 peptide to the TAR RNA target. The OPES simulation delineated the mechanism of peptide binding and unbinding to and from the RNA and enabled the calculation of the underlying free energy landscape. Our results demonstrate the potential of the Deep-TDA method for designing CVs to study complex biomolecular recognition processes.

Enhanced Sampling with Suboptimal Collective Variables: Reconciling Accuracy and Convergence Speed

We introduce an enhanced sampling algorithm to obtain converged free energy landscapes of molecular rare events, even when the collective variable (CV) used for biasing is not optimal. Our approach samples a time-dependent target distribution by combining the on-the-fly probability enhanced sampling and its exploratory variant, OPES Explore. This promotes more transitions between the relevant metastable states and accelerates the convergence speed of the free energy estimate. We demonstrate the successful application of this combined algorithm on the two-dimensional Wolfe-Quapp potential, millisecond time-scale ligand-receptor binding in the trypsin-benzamidine complex, and folding-unfolding transitions in chignolin mini-protein. Our proposed algorithm can compute accurate free energies at an affordable computational cost and is robust in terms of the choice of CVs, making it particularly promising for the simulation of complex biomolecular systems.

Descriptor-Free Collective Variables from Geometric Graph Neural Networks

Enhanced sampling simulations make the computational study of rare events feasible. A large family of such methods crucially depends on the definition of some collective variables (CVs) that could provide a low-dimensional representation of the relevant physics of the process. Recently, many methods have been proposed to semiautomatize the CV design by using machine learning tools to learn the variables directly from the simulation data. However, most methods are based on feedforward neural networks and require some user-defined physical descriptors. Here, we propose bypassing this step using a graph neural network to directly use the atomic coordinates as input for the CV model. This way, we achieve a fully automatic approach to CV determination that provides variables invariant under the relevant symmetries, especially the permutational one. Furthermore, we provide different analysis tools to favor the physical interpretation of the final CV. We prove the robustness of our approach using different methods from the literature for the optimization of the CV, and we prove its efficacy on several systems, including a small peptide, an ion dissociation in explicit solvent, and a simple chemical reaction.

Methane degassing in global river reservoirs and its impacts on carbon budgets and sustainable water management

Degassing methane (CH4) through reservoir water compromises hydroelectricity's presumed low-carbon status, which has emerged as a critical hotspot for global carbon dynamics. However, a comprehensive understanding of the involved pathways remains elusive, hindering the accurate estimation of global reservoirs' carbon budget (emission-to-burial ratio). This study presents a holistic upscaling approach to assess methane degassing in global river reservoirs and its impacts on carbon budgets. Firstly, a machine learning model is utilized to characterize the contributions of climate and human factors to annual water residence time. Secondly, the stepwise multiple linear regression method is used to calculate CH4 degassing emissions for each reservoir. Finally, to systematically tackle all sources of uncertainty, separate uncertainty analyses are implemented for the estimates of degassing emissions, areal emissions, and organic carbon burial. Analyzing 30-year data from 6695 reservoirs worldwide, our assessment considers water residence time, temperature, catchment area, and reservoir size. Findings indicate that water releases contribute significantly to global CO2 emissions from reservoirs, elevating the carbon budget by 20 % from 2.02 to 2.18 TgC/year, underscoring the previously underestimated significance of CH4 degassing in shaping the carbon cycle impact of river reservoirs. We propose a redefined threshold for low carbon credits, suggesting that reservoirs with power densities exceeding 6.1 MW/km2, instead of the conventional 4 MW/km2, should qualify. This study underscores the need for sustainable water management and reshaping the carbon dynamics associated with hydroelectricity. Future research can advocate Artificial Intelligence (AI) techniques to enhance water management and mitigate carbon emissions by multi-objectively optimizing reservoir operations.

Correlating enzymatic reactivity for different substrates using transferable data-driven collective variables

Machine learning (ML) is transforming the investigation of complex biological processes. In enzymatic catalysis, one significant challenge is identifying the reactive conformations (RC) of the enzyme:substrate complex where the substrate assumes a precise arrangement in the active site necessary to initiate a reaction. Traditional methods are hindered by the complexity of the multidimensional free energy landscape involved in the transition from nonreactive to reactive conformations. Here, we applied ML techniques to address this challenge, focusing on human pancreatic α-amylase, a crucial enzyme in type-II diabetes treatment. Using ML-based collective variables (CVs), we correlated the probability of being in a RC with the experimental catalytic activity of several malto-oligosaccharide substrates. Our findings demonstrate a remarkable transferability of these CVs across various compounds, significantly streamlining the modeling process and reducing both computational demand and manual intervention in setting up simulations for new substrates. This approach not only advances our understanding of enzymatic processes but also holds substantial potential for accelerating drug discovery by enabling rapid and accurate evaluation of drug efficacy across different generations of inhibitors.

Preparation of micron-sized benzamidine-modified magnetic agarose beads for trypsin purification from fish viscera

The effective method for trypsin purification should be established because trypsin has important economic value. In this work, a novel and simple strategy was proposed for fabricating micron-sized magnetic Fe3O4@agarose-benzamidine beads (MABB) with benzamidine as a ligand, which can efficiently and selectively capture trypsin. The micro-sized MABB, with clear spherical core-shell structure and average particle size of 6.6 μm, showed excellent suspension ability and magnetic responsiveness in aqueous solution. The adsorption capacity and selectivity of MABB towards target trypsin were significantly better than those of non-target lysozyme. According to the Langmuir equation, the maximum adsorption capacity of MABB for trypsin was 1946 mg g-1 at 25 °C, and the adsorption should be a physical sorption process. Furthermore, the initial adsorption rate and half equilibrium time of MABB toward trypsin were 787.4 mg g-1 min-1 and 0.71 min, respectively. To prove the practicability, MABB-based magnetic solid-phase extraction (MSPE) was proposed, and the related parameters were optimized in detail to improve the purification efficiency. With Tris-HCl buffer (50 mM, 10 mM CaCl2, pH 8.0) as extraction buffer, Tris-HCl buffer (50 mM, 100 mM CaCl2, pH 8.0) as rinsing buffer, acidic eluent (0.01 M HCl, 0.5 M NaCl, pH 2.0) as eluent buffer and alkaline buffer (1 M Tris-HCl buffer, pH 10.0) as neutralization solution, the MABB-based MSPE was successfully used for trypsin purification from the viscera of grass carp (Ctenopharyngodon idella). The molecular weight of purified trypsin was determined as approximate 23 kDa through sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE). The purified trypsin was highly active from 30 °C to 60 °C, with an optimum temperature of 50 °C, and was tolerant to pH variation, exhibiting 85 % of maximum enzyme activity from pH 7.0 to 10.0. The results demonstrated that the proposed MABB-based MSPE could effectively purify trypsin and ensure the biological activity of purified trypsin. Therefore, we believe that the novel MABB could be applicable for efficient purification of trypsin from complex biological systems.

Advances and Challenges in Milestoning Simulations for Drug-Target Kinetics

Molecular dynamics simulations have become indispensable for exploring complex biological processes, yet their limitations in capturing rare events hinder our understanding of drug-target kinetics. In this Perspective, we investigate the domain of milestoning simulations to understand this challenge. The milestoning approach divides the phase space of the drug-target complex into discrete cells, offering extended time scale insights. This Perspective traces the history, applications, and future potential of milestoning simulations in the context of drug-target kinetics. It explores the fundamental principles of milestoning, highlighting the importance of probabilistic transitions and transition time independence. Markovian milestoning with Voronoi tessellations is revisited to address the traditional milestoning challenges. While observing the advancements in this field, this Perspective also addresses impending challenges in estimating drug-target unbinding rate constants through milestoning simulations, paving the way for more effective drug design strategies.

Host-Guest Binding Free Energies à la Carte: An Automated OneOPES Protocol

Estimating absolute binding free energies from molecular simulations is a key step in computer-aided drug design pipelines, but the agreement between computational results and experiments is still very inconsistent. Both the accuracy of the computational model and the quality of the statistical sampling contribute to this discrepancy, yet disentangling the two remains a challenge. In this study, we present an automated protocol based on OneOPES, an enhanced sampling method that exploits replica exchange and can accelerate several collective variables to address the sampling problem. We apply this protocol to 37 host-guest systems. The simplicity of setting up the simulations and producing well-converged binding free energy estimates without the need to optimize simulation parameters provides a reliable solution to the sampling problem. This, in turn, allows for a systematic force field comparison and ranking according to the correlation between simulations and experiments, which can inform the selection of an appropriate model. The protocol can be readily adapted to test more force field combinations and study more complex protein-ligand systems, where the choice of an appropriate physical model is often based on heuristic considerations rather than systematic optimization.

Absolute Binding Free Energies with OneOPES

The calculation of absolute binding free energies (ABFEs) for protein-ligand systems has long been a challenge. Recently, refined force fields and algorithms have improved the quality of the ABFE calculations. However, achieving the level of accuracy required to inform drug discovery efforts remains difficult. Here, we present a transferable enhanced sampling strategy to accurately calculate absolute binding free energies using OneOPES with simple geometric collective variables. We tested the strategy on two protein targets, BRD4 and Hsp90, complexed with a total of 17 chemically diverse ligands, including both molecular fragments and drug-like molecules. Our results show that OneOPES accurately predicts protein-ligand binding affinities with a mean unsigned error within 1 kcal mol-1 of experimentally determined free energies, without the need to tailor the collective variables to each system. Furthermore, our strategy effectively samples different ligand binding modes and consistently matches the experimentally determined structures regardless of the initial protein-ligand configuration. Our results suggest that the proposed OneOPES strategy can be used to inform lead optimization campaigns in drug discovery and to study protein-ligand binding and unbinding mechanisms.

Bidentate Coordination Structure Facilitates High-Voltage and High-Utilization Aqueous Zn-I2 Batteries

The aqueous zinc-iodine battery is a promising energy storage device, but the conventional two-electron reaction potential and energy density of the iodine cathode are far from meeting practical application requirements. Given that iodine is rich in redox reactions, activating the high-valence iodine cathode reaction has become a promising research direction for developing high-voltage zinc-iodine batteries. In this work, by designing a multifunctional electrolyte additive trimethylamine hydrochloride (TAH), a stable high-valence iodine cathode in four-electron-transfer I-/I2/I+ reactions with a high theoretical specific capacity is achieved through a unique amine group, Cl bidentate coordination structure of (TA)ICl. Characterization techniques such as synchrotron radiation, in situ Raman spectra, and DFT calculations are used to verify the mechanism of the stable bidentate structure. This electrolyte additive stabilizes the zinc anode by promoting the desolvation process and shielding mechanism, enabling the zinc anode to cycle steadily at a maximum areal capacity of 57 mAh cm-2 with 97 % zinc utilization rate. Finally, the four-electron-transfer aqueous Zn-I2 full cell achieves 5000 stable cycles at an N/P ratio of 2.5. The unique bidentate coordination structure contributes to the further development of high-valence and high capacity aqueous zinc-iodine batteries.

Incorporating Prior Knowledge in the Seeds of Adaptive Sampling Molecular Dynamics Simulations of Ligand Transport in Enzymes with Buried Active Sites

Because most proteins have buried active sites, protein tunnels or channels play a crucial role in the transport of small molecules into buried cavities for enzymatic catalysis. Tunnels can critically modulate the biological process of protein-ligand recognition. Various molecular dynamics methods have been developed for exploring and exploiting the protein-ligand conformational space to extract high-resolution details of the binding processes, a recent example being energetically unbiased high-throughput adaptive sampling simulations. The current study systematically contrasted the role of integrating prior knowledge while generating useful initial protein-ligand configurations, called seeds, for these simulations. Using a nontrivial system of a haloalkane dehalogenase mutant with multiple transport tunnels leading to a deeply buried active site, simulations were employed to derive kinetic models describing the process of association and dissociation of the substrate molecule. The most knowledge-based seed generation enabled high-throughput simulations that could more consistently capture the entire transport process, explore the complex network of transport tunnels, and predict equilibrium dissociation constants, koff/kon, on the same order of magnitude as experimental measurements. Overall, the infusion of more knowledge into the initial seeds of adaptive sampling simulations could render analyses of transport mechanisms in enzymes more consistent even for very complex biomolecular systems, thereby promoting drug development efforts and the rational design of enzymes with buried active sites.

Do All Paths Lead to Rome? How Reliable is Umbrella Sampling Along a Single Path?

Molecular dynamics (MD) simulations are widely applied to estimate absolute binding free energies of protein-ligand and protein-protein complexes. A routinely used method for binding free energy calculations with MD is umbrella sampling (US), which calculates the potential of mean force (PMF) along a single reaction coordinate. Surprisingly, in spite of its widespread use, few validation studies have focused on the convergence of the free energy computed along a single path for specific cases, not addressing the reproducibility of such calculations in general. In this work, we therefore investigate the reproducibility and convergence of US along a standard distance-based reaction coordinate for various protein-protein and protein-ligand complexes, following commonly used guidelines for the setup. We show that repeating the complete US workflow can lead to differences of 2-20 kcal/mol in computed binding free energies. We attribute those discrepancies to small differences in the binding pathways. While these differences are unavoidable in the established US protocol, the popularity of the latter could hint at a lack of awareness of such reproducibility problems. To test if the convergence of PMF profiles can be improved if multiple pathways are sampled simultaneously, we performed additional simulations with an adaptive-biasing method, here the accelerated weight histogram (AWH) approach. Indeed, the PMFs obtained from AHW simulations are consistent and reproducible for the systems tested. To the best of our knowledge, our work is the first to attempt a systematic assessment of the pitfalls in one the most widely used protocols for computing binding affinities. We anticipate therefore that our results will provide an incentive for a critical reassessment of the validity of PMFs computed with US, and make a strong case to further benchmark the performance of adaptive-biasing methods for computing binding affinities.

Machine Learning Derived Collective Variables for the Study of Protein Homodimerization in Membrane

The accurate calculation of equilibrium constants for protein-protein association is of fundamental importance to quantitative biology and remains an outstanding challenge for computational biophysics. Traditionally, equilibrium constants have been computed from one-dimensional free energy surfaces derived from sampling along a single collective variable. Importantly, recent advances in enhanced sampling methodology have facilitated the characterization of multidimensional free energy landscapes, often exposing multiple thermodynamically important minima missed by more restrictive sampling methods. A key to the effectiveness of this multidimensional sampling approach is the identification of collective variables that effectively define the configurational space of dissociated and associated states. Here we present the application of two machine learning methods for the unbiased determination of collective variables for enhanced sampling and analysis of protein-protein association. Our results both validate prior work, based on intuition derived collective variables, and demonstrate the effectiveness of the machine learning methods for the identification of collective variables for association reactions in complex biomolecular systems.

Leveraging Machine Learning-Guided Molecular Simulations Coupled with Experimental Data to Decipher Membrane Binding Mechanisms of Aminosterols

Understanding the molecular mechanisms of the interactions between specific compounds and cellular membranes is essential for numerous biotechnological applications, including targeted drug delivery, elucidation of the drug mechanism of action, pathogen identification, and novel antibiotic development. However, estimation of the free energy landscape associated with solute binding to realistic biological systems is still a challenging task. In this work, we leverage the Time-lagged Independent Component Analysis (TICA) in combination with neural networks (NN) through the Deep-TICA approach for determining the free energy associated with the membrane insertion processes of two natural aminosterol compounds, trodusquemine (TRO), and squalamine (SQ). These compounds are particularly noteworthy because they interact with the outer layer of neuron membranes, protecting them from the toxic action of misfolded proteins involved in neurodegenerative disorders, in both their monomeric and oligomeric forms. We demonstrate how this strategy could be used to generate an effective collective variable for describing solute absorption in the membrane and for estimating free energy landscape of translocation via on-the-fly probability enhanced sampling (OPES) method. In this context, the computational protocol allowed an exhaustive characterization of the aminosterol entry pathway into a neuron-like lipid bilayer. Furthermore, it provided accurate prediction of membrane binding affinities, in close agreement with the experimental binding data obtained by using fluorescently labeled aminosterols and large unilamellar vesicles (LUVs). The findings contribute significantly to our understanding of aminosterol entry pathways and aminosterol-lipid membrane interactions. Finally, the computational methods deployed in this study further demonstrate considerable potential for investigating membrane binding processes.

Rotational Dynamics of Water near Osmolytes by Molecular Dynamics Simulations

The behavior of water molecules around organic molecules has attracted considerable attention as a crucial factor influencing the properties and functions of soft matter and biomolecules. Recently, it has been suggested that the change in protein stability upon the addition of small organic molecules (osmolytes) is dominated by the change in the water dynamics caused by the osmolyte, where the dynamics of not only the directly interacting water molecules but also the long-range hydration layer affect the protein stability. However, the relation between the long-range structure of hydration water in various solutions and the water dynamics remains unclear at the molecular level. We performed density-functional tight-binding molecular dynamics simulations to elucidate the varying rotational dynamics of water molecules in 15 osmolyte solutions. A positive correlation was observed between the rotational relaxation time and our proposed normalized parameter obtained by dividing the number of hydrogen bonds between water molecules by the number of nearest-neighbor water molecules. For the 15 osmolyte solutions, an increase or a decrease in the value of the normalized parameter for the second hydration shell tended to result in water molecules with slow and fast rotational dynamics, respectively, thus illustrating the importance of the second hydration shell for the rotational dynamics of water molecules. Our simulation results are anticipated to advance the current understanding of water dynamics around organic molecules and the long-range structure of water molecules.

Data-driven classification of ligand unbinding pathways

Studying the pathways of ligand-receptor binding is essential to understand the mechanism of target recognition by small molecules. The binding free energy and kinetics of protein-ligand complexes can be computed using molecular dynamics (MD) simulations, often in quantitative agreement with experiments. However, only a qualitative picture of the ligand binding/unbinding paths can be obtained through a conventional analysis of the MD trajectories. Besides, the higher degree of manual effort involved in analyzing pathways limits its applicability in large-scale drug discovery. Here, we address this limitation by introducing an automated approach for analyzing molecular transition paths with a particular focus on protein-ligand dissociation. Our method is based on the dynamic time-warping algorithm, originally designed for speech recognition. We accurately classified molecular trajectories using a very generic descriptor set of contacts or distances. Our approach outperforms manual classification by distinguishing between parallel dissociation channels, within the pathways identified by visual inspection. Most notably, we could compute exit-path-specific ligand-dissociation kinetics. The unbinding timescale along the fastest path agrees with the experimental residence time, providing a physical interpretation to our entirely data-driven protocol. In combination with appropriate enhanced sampling algorithms, this technique can be used for the initial exploration of ligand-dissociation pathways as well as for calculating path-specific thermodynamic and kinetic properties.

Thermodynamically Optimized Machine-Learned Reaction Coordinates for Hydrophobic Ligand Dissociation

Ligand unbinding is mediated by its free energy change, which has intertwined contributions from both energy and entropy. It is important, but not easy, to quantify their individual contributions to the free energy profile. We model hydrophobic ligand unbinding for two systems, a methane particle and a C60 fullerene, both unbinding from hydrophobic pockets in all-atom water. Using a modified deep learning framework, we learn a thermodynamically optimized reaction coordinate to describe the hydrophobic ligand dissociation for both systems. Interpretation of these reaction coordinates reveals the roles of entropic and enthalpic forces as the ligand and pocket sizes change. In both cases, we observe that the free-energy barrier to unbinding is dominated by entropy considerations. Furthermore, the process of methane unbinding is driven by methane solvation, while fullerene unbinding is driven first by pocket wetting and then fullerene wetting. For both solutes, the direct importance of the distance from the binding pocket to the learned reaction coordinate is present, but low. Our framework and subsequent feature important analysis thus give useful thermodynamic insight into hydrophobic ligand dissociation problems that are otherwise difficult to glean.

Cation Chloride Cotransporter NKCC1 Operates through a Rocking-Bundle Mechanism

The sodium, potassium, and chloride cotransporter 1 (NKCC1) plays a key role in tightly regulating ion shuttling across cell membranes. Lately, its aberrant expression and function have been linked to numerous neurological disorders and cancers, making it a novel and highly promising pharmacological target for therapeutic interventions. A better understanding of how NKCC1 dynamically operates would therefore have broad implications for ongoing efforts toward its exploitation as a therapeutic target through its modulation. Based on recent structural data on NKCC1, we reveal conformational motions that are key to its function. Using extensive deep-learning-guided atomistic simulations of NKCC1 models embedded into the membrane, we captured complex dynamical transitions between alternate open conformations of the inner and outer vestibules of the cotransporter and demonstrated that NKCC1 has water-permeable states. We found that these previously undefined conformational transitions occur via a rocking-bundle mechanism characterized by the cooperative angular motion of transmembrane helices (TM) 4 and 9, with the contribution of the extracellular tip of TM 10. We found these motions to be critical in modulating ion transportation and in regulating NKCC1's water transporting capabilities. Specifically, we identified interhelical dynamical contacts between TM 10 and TM 6, which we functionally validated through mutagenesis experiments of 4 new targeted NKCC1 mutants. We conclude showing that those 4 residues are highly conserved in most Na+-dependent cation chloride cotransporters (CCCs), which highlights their critical mechanistic implications, opening the way to new strategies for NKCC1's function modulation and thus to potential drug action on selected CCCs.

Roles of Accelerated Molecular Dynamics Simulations in Predictions of Binding Kinetic Parameters

Rational predictions on binding kinetics parameters of drugs to targets play significant roles in future drug designs. Full conformational samplings of targets are requisite for accurate predictions of binding kinetic parameters. In this review, we mainly focus on the applications of enhanced sampling technologies in calculations of binding kinetics parameters and residence time of drugs. The methods involved in molecular dynamics simulations are applied to not only probe conformational changes of targets but also reveal calculations of residence time that is significant for drug efficiency. For this review, special attention are paid to accelerated molecular dynamics (aMD) and Gaussian aMD (GaMD) simulations that have been adopted to predict the association or disassociation rate constant. We also expect that this review can provide useful information for future drug design.

Perspectives on Ligand/Protein Binding Kinetics Simulations: Force Fields, Machine Learning, Sampling, and User-Friendliness

Computational techniques applied to drug discovery have gained considerable popularity for their ability to filter potentially active drugs from inactive ones, reducing the time scale and costs of preclinical investigations. The main focus of these studies has historically been the search for compounds endowed with high affinity for a specific molecular target to ensure the formation of stable and long-lasting complexes. Recent evidence has also correlated the in vivo drug efficacy with its binding kinetics, thus opening new fascinating scenarios for ligand/protein binding kinetic simulations in drug discovery. The present article examines the state of the art in the field, providing a brief summary of the most popular and advanced ligand/protein binding kinetics techniques and evaluating their current limitations and the potential solutions to reach more accurate kinetic models. Particular emphasis is put on the need for a paradigm change in the present methodologies toward ligand and protein parametrization, the force field problem, characterization of the transition states, the sampling issue, and algorithms' performance, user-friendliness, and data openness.

Kinetics from Metadynamics: Principles, Applications, and Outlook

Metadynamics is a popular enhanced sampling algorithm for computing the free energy landscape of rare events by using molecular dynamics simulation. Ten years ago, Tiwary and Parrinello introduced the infrequent metadynamics approach for calculating the kinetics of transitions across free energy barriers. Since then, metadynamics-based methods for obtaining rate constants have attracted significant attention in computational molecular science. Such methods have been applied to study a wide range of problems, including protein-ligand binding, protein folding, conformational transitions, chemical reactions, catalysis, and nucleation. Here, we review the principles of elucidating kinetics from metadynamics-like approaches, subsequent methodological developments in this area, and successful applications on chemical, biological, and material systems. We also highlight the challenges of reconstructing accurate kinetics from enhanced sampling simulations and the scope of future developments.

OneOPES, a Combined Enhanced Sampling Method to Rule Them All

Enhanced sampling techniques have revolutionized molecular dynamics (MD) simulations, enabling the study of rare events and the calculation of free energy differences in complex systems. One of the main families of enhanced sampling techniques uses physical degrees of freedom called collective variables (CVs) to accelerate a system's dynamics and recover the original system's statistics. However, encoding all the relevant degrees of freedom in a limited number of CVs is challenging, particularly in large biophysical systems. Another category of techniques, such as parallel tempering, simulates multiple replicas of the system in parallel, without requiring CVs. However, these methods may explore less relevant high-energy portions of the phase space and become computationally expensive for large systems. To overcome the limitations of both approaches, we propose a replica exchange method called OneOPES that combines the power of multireplica simulations and CV-based enhanced sampling. This method efficiently accelerates the phase space sampling without the need for ideal CVs, extensive parameters fine tuning nor the use of a large number of replicas, as demonstrated by its successful applications to protein-ligand binding and protein folding benchmark systems. Our approach shows promise as a new direction in the development of enhanced sampling techniques for molecular dynamics simulations, providing an efficient and robust framework for the study of complex and unexplored problems.

Investigating the Unbinding of Muscarinic Antagonists from the Muscarinic 3 Receptor

Patient symptom relief is often heavily influenced by the residence time of the inhibitor-target complex. For the human muscarinic receptor 3 (hMR3), tiotropium is a long-acting bronchodilator used in conditions such as asthma or chronic obstructive pulmonary disease (COPD). The mechanistic insights into this inhibitor remain unclear; specifically, the elucidation of the main factors determining the unbinding rates could help develop the next generation of antimuscarinic agents. Using our novel unbinding algorithm, we were able to investigate ligand dissociation from hMR3. The unbinding paths of tiotropium and two of its analogues, N-methylscopolamin and homatropine methylbromide, show a consistent qualitative mechanism and allow us to identify the structural bottleneck of the process. Furthermore, our machine learning-based analysis identified key roles of the ECL2/TM5 junction involved in the transition state. Additionally, our results point to relevant changes at the intracellular end of the TM6 helix leading to the ICL3 kinase domain, highlighting the closest residue L482. This residue is located right between two main protein binding sites involved in signal transduction for hMR3's activation and regulation. We also highlight key pharmacophores of tiotropium that play determining roles in the unbinding kinetics and could aid toward drug design and lead optimization.

Improving binding entropy by higher ligand symmetry? - A case study with human matriptase

Understanding different contributions to the binding entropy of ligands is of utmost interest to better predict affinity and the thermodynamic binding profiles of protein-ligand interactions and to develop new strategies for ligand optimization. To these means, the largely neglected effects of introducing higher ligand symmetry, thereby reducing the number of energetically distinguishable binding modes on binding entropy using the human matriptase as a model system, were investigated. A set of new trivalent phloroglucinol-based inhibitors that address the roughly symmetric binding site of the enzyme was designed, synthesized, and subjected to isothermal titration calorimetry. These highly symmetric ligands that can adopt multiple indistinguishable binding modes exhibited high entropy-driven affinity in line with affinity-change predictions.

Predicting Protein-Ligand Binding and Unbinding Kinetics with Biased MD Simulations and Coarse-Graining of Dynamics: Current State and Challenges

The prediction of drug-target binding and unbinding kinetics that occur on time scales between milliseconds and several hours is a prime challenge for biased molecular dynamics simulation approaches. This Perspective gives a concise summary of the theory and the current state-of-the-art of such predictions via biased simulations, of insights into the molecular mechanisms defining binding and unbinding kinetics as well as of the extraordinary challenges predictions of ligand kinetics pose in comparison to binding free energy predictions.

Predicting Biomolecular Binding Kinetics: A Review

Biomolecular binding kinetics including the association (k_on) and dissociation (k_off) rates are critical parameters for therapeutic design of small-molecule drugs, peptides, and antibodies. Notably, the drug molecule residence time or dissociation rate has been shown to correlate with their efficacies better than binding affinities. A wide range of modeling approaches including quantitative structure-kinetic relationship models, Molecular Dynamics simulations, enhanced sampling, and Machine Learning has been developed to explore biomolecular binding and dissociation mechanisms and predict binding kinetic rates. Here, we review recent advances in computational modeling of biomolecular binding kinetics, with an outlook for future improvements.

Advances in computational methods for ligand binding kinetics

Binding kinetic parameters can be correlated with drug efficacy, which in recent years led to the development of various computational methods for predicting binding kinetic rates and gaining insight into protein-drug binding paths and mechanisms. In this review, we introduce and compare computational methods recently developed and applied to two systems, trypsin-benzamidine and kinase-inhibitor complexes. Methods involving enhanced sampling in molecular dynamics simulations or machine learning can be used not only to predict kinetic rates, but also to reveal factors modulating the duration of residence times, selectivity, and drug resistance to mutations. Methods which require less computational time to make predictions are highlighted, and suggestions to reduce the error of computed kinetic rates are presented.

Kinetics of Drug Release from Clay Using Enhanced Sampling Methods

A key step in the development of a new drug, is the design of drug-excipient complexes that lead to optimal drug release kinetics. Computational chemistry and specifically enhanced sampling molecular dynamics methods can play a key role in this context, by minimizing the need for expensive experiments, and reducing cost and time. Here we show that recent advances in enhanced sampling methodologies can be brought to fruition in this area. We demonstrate the potential of these methodologies by simulating the drug release kinetics of the complex praziquantel-montmorillonite in water. Praziquantel finds promising applications in the treatment of schistosomiasis, but its biopharmaceutical profile needs to be improved, and a cheap material such as the montmorillonite clay would be a very convenient excipient. We simulate the drug release both from surface and interlayer space, and find that the diffusion of the praziquantel inside the interlayer space is the process that limits the rate of drug release.

Rare Event Kinetics from Adaptive Bias Enhanced Sampling

We introduce a novel enhanced sampling approach named on-the-fly probability enhanced sampling (OPES) flooding for calculating the kinetics of rare events from atomistic molecular dynamics simulation. This method is derived from the OPES approach [Invernizzi and Parrinello, J. Phys. Chem. Lett. 2020, 11, 7, 2731-2736], which has been recently developed for calculating converged free energy surfaces for complex systems. In this paper, we describe the theoretical details of the OPES flooding technique and demonstrate the application on three systems of increasing complexity: barrier crossing in a two-dimensional double-well potential, conformational transition in the alanine dipeptide in the gas phase, and the folding and unfolding of the chignolin polypeptide in an aqueous environment. From extensive tests, we show that the calculation of accurate kinetics not only requires the transition state to be bias-free, but the amount of bias deposited should also not exceed the effective barrier height measured along the chosen collective variables. In this vein, the possibility of computing rates from biasing suboptimal order parameters has also been explored. Furthermore, we describe the choice of optimum parameter combinations for obtaining accurate results from limited computational effort.