Taming Rugged Free Energy Landscapes Using an Average Force

Complementary Peptide Interactions Support the Ultra-Rigidity of Polymers of De Novo Designed Click-Functionalized Bundlemers

Computationally designed 29-residue peptides yield tetra-α-helical bundles with D2 symmetry. The "bundlemers" can be bifunctionally linked via thiol-maleimide cross-links at their N-termini, yielding supramolecular polymers with unusually large, micrometer-scale persistence lengths. To provide a molecularly resolved understanding of these systems, all-atom molecular modeling and simulations of linked bundlemers in explicit solvent are presented. A search over relative orientations of the bundlemers identifies a structure, wherein at the bundlemer-bundlemer interface, interior hydrophobic residues are in contact, and α-helices are aligned with a pseudocontiguous α-helix that spans the interface. Calculation of a potential of mean force confirms that the structure in which the bundlemers are in contact and colinearly aligned is a stable minimum. Analyses of hydrogen bonds and hydrophobic complementarity highlight the complementary interactions at the interface. The molecular insight provided reveals the molecular origins of bundlemer alignment within the supramolecular polymers.

Probing the orientation and membrane permeation of rhodamine voltage reporters through molecular simulations and free energy calculations

The transmembrane potential of plasma membranes and membrane-bound organelles plays a fundamental role in cellular functions such as signal transduction, ATP synthesis, and homeostasis. Rhodamine voltage reporters (RhoVRs), which operate based on the photoinduced electron transfer (PeT) mechanism, are non-invasive, small-molecule voltage sensors that can detect rapid voltage changes, with some of them specifically targeting the inner mitochondrial membrane. In this work, we conducted extensive molecular dynamics simulations and free-energy calculations to investigate the physicochemical properties governing the orientation as well as membrane permeation barriers of three RhoVRs. Our results indicate that the positioning of the most polarized functional group relative to the hydrophobic molecular wire dictates the alignment of RhoVRs with the membrane normal, thereby, significantly affecting their voltage sensitivity. Free-energy calculations in different membrane systems identify significantly higher barriers against the permeation of RhoVR 1 compared to SPIRIT RhoVR 1, explaining their distinct subcellular localization profiles. Subsequent free-energy calculations of the distinguishing components from the two different RhoVRs provide additional insight into the physicochemical properties governing their membrane permeation. The connection between chemical composition and membrane orientation, as well as permeation behaviors of RhoVRs revealed by our calculations provides general guiding principles for the rational design of PeT-based fluorescent dyes with enhanced voltage sensitivity and desired subcellular distribution.

Expanded Functionality and Portability for the Colvars Library

Colvars is an open-source C++ library that provides a modular toolkit for collective-variable-based molecular simulations. It allows practitioners to easily create and implement descriptors that best fit a process of interest and to apply a wide range of biasing algorithms in collective variable space. This paper reviews several features and improvements to Colvars that were added since its original introduction. Special attention is given to contributions that significantly expanded the capabilities of this software or its distribution with major MD simulation packages. Collective variables can now be optimized either manually or by machine-learning methods, and the space of descriptors can be explored interactively using the graphical interface included in VMD. Beyond the spatial coordinates of individual molecules, Colvars can now apply biasing forces to mesoscale structures and alchemical degrees of freedom and perform simulations guided by experimental data within ensemble averages or probability distributions. It also features advanced computational schemes to boost the accuracy, robustness, and general applicability of simulation methods, including extended-system and multiple-walker adaptive biasing force, boundary conditions for metadynamics, replica exchange with biasing potentials, and adiabatic bias molecular dynamics. The library is made available directly within the main distributions of the academic software GROMACS, LAMMPS, NAMD, Tinker-HP, and VMD. The robustness of the software and the reliability of the results are ensured through the use of continuous integration with a test suite within the source repository.

Isoleucine gate blocks K+ conduction in C-type inactivation

Many voltage-gated potassium (Kv) channels display a time-dependent phenomenon called C-type inactivation, whereby prolonged activation by voltage leads to the inhibition of ionic conduction, a process that involves a conformational change at the selectivity filter toward a non-conductive state. Recently, a high-resolution structure of a strongly inactivated triple-mutant channel kv1.2-kv2.1-3m revealed a novel conformation of the selectivity filter that is dilated at its outer end, distinct from the well-characterized conductive state. While the experimental structure was interpreted as the elusive non-conductive state, our molecular dynamics simulations and electrophysiological measurements show that the dilated filter of kv1.2-kv2.1-3m is conductive and, as such, cannot completely account for the inactivation of the channel observed in the structural experiments. The simulation shows that an additional conformational change, implicating isoleucine residues at position 398 along the pore lining segment S6, is required to effectively block ion conduction. The I398 residues from the four subunits act as a state-dependent hydrophobic gate located immediately beneath the selectivity filter. These observations are corroborated by electrophysiological experiments showing that ion permeation can be resumed in the kv1.2-kv2.1-3m channel when I398 is mutated to an asparagine-a mutation that does not abolish C-type inactivation since digitoxin (AgTxII) fails to block the ionic permeation of kv1.2-kv2.1-3m_I398N. As a critical piece of the C-type inactivation machinery, this structural feature is the potential target of a broad class of quaternary ammonium (QA) blockers and negatively charged activators thus opening new research directions toward the development of drugs that specifically modulate gating states of Kv channels.

Cobaltabis(dicarbollide) Interaction with DNA Resolved at the Atomic Scale

Boron neutron capture therapy represents a promising avenue for cancer treatment that requires nontoxic drugs with a high boron content efficiently distributed into cancerous cells. The metallacarborane o-cobaltabis(dicarbollide) ([COSAN]-) fulfills these requirements and constitutes an attractive candidate. Nevertheless, the interaction of this promising drug with nucleic acids, the assumed target of the biological damage, is poorly understood since contradictory results are reported in the literature. This work establishes the DNA/[COSAN]- interaction strength, mechanism, and time scale at the atomistic level by using a combination of microsecond-molecular dynamics and hybrid quantum mechanics/molecular mechanics simulations and by quantifying the absolute binding free energy. Results show that the DNA/[COSAN]- interaction is highly dependent on the ionic strength of the medium. A relatively weak DNA major groove binding (ΔGbind= -2.49 kcal/mol) driven mostly by dihydrogen B-H···H-N bonding is observed in the simulations only at a high NaCl concentration, whereas DNA intercalation mode is deemed highly unlikely.

Dissected antiporter modules establish minimal proton-conduction elements of the respiratory complex I

The respiratory Complex I is a highly intricate redox-driven proton pump that powers oxidative phosphorylation across all domains of life. Yet, despite major efforts in recent decades, its long-range energy transduction principles remain highly debated. We create here minimal proton-conducting membrane modules by engineering and dissecting the key elements of the bacterial Complex I. By combining biophysical, biochemical, and computational experiments, we show that the isolated antiporter-like modules of Complex I comprise all functional elements required for conducting protons across proteoliposome membranes. We find that the rate of proton conduction is controlled by conformational changes of buried ion-pairs that modulate the reaction barriers by electric field effects. The proton conduction is also modulated by bulky residues along the proton channels that are key for establishing a tightly coupled proton pumping machinery in Complex I. Our findings provide direct experimental evidence that the individual antiporter modules are responsible for the proton transport activity of Complex I. On a general level, our findings highlight electrostatic and conformational coupling mechanisms in the modular energy-transduction machinery of Complex I with distinct similarities to other enzymes.

Mechanism of Dual-Site Recognition in a Classic DNA Aptamer

Nucleic acid aptamers possess unique advantages in specific recognition. However, the lack of in-depth investigation into their dynamic recognition mechanisms has restricted their rational design and potential applications in fields such as biosensing and targeted therapy. We herein utilized enhanced sampling molecular dynamics to address affinities of adenosine monophosphate (AMP) to the dual binding sites in the DNA aptamer, focusing on the dynamic recognition mechanism and pathways. The present results indicate that in addition to the widely known intermolecular interactions, inequivalence of chemical environments of the two binding sites leads to slightly higher stability of AMP binding to the site proximal to the aptamer terminus. In the presence of two AMPs captured by the two sites, each binding free energy is enhanced. In particular, an additional hydrogen bond of AMP to A10 is introduced in the dual-site binding complex, which increases the binding energy from -4.25 ± 0.47 to -9.48 ± 0.33 kcal mol-1 in the site close to the loop. For the dual-site recognition process, the free energy landscape and minimum free energy pathway calculations elucidate the crucial role of electrostatic interactions between the AMP phosphate groups and Na+ ions in positively cooperative binding mechanisms.

Overcoming Sampling Issues and Improving Computational Efficiency in Collective-Variable-Based Enhanced-Sampling Simulations: A Tutorial

This tutorial is designed to help users overcome sampling challenges and improve computational efficiency in collective-variable (CV)-based enhanced-sampling, or importance-sampling, simulations. Toward this end, we introduce well-tempered metadynamics-extended adaptive biasing force (WTM-eABF) and its integration with Gaussian accelerated molecular dynamics (GaMD). Additionally, use will be made of a method for identifying the least-free-energy pathway (LFEP) and multiple concurrent pathways on high-dimensional free-energy surfaces. We illustrate these sampling techniques with the conformational equilibria of trialanine and chignolin in aqueous solution as test cases. This tutorial assumes that the user has prior experience with molecular dynamics (MD) simulations, in general, with the popular program NAMD, and to some extent with Colvars, the module for CV-based calculations. This tutorial can, however, in large measure be used in conjunction with alternate MD engines that support the Colvars module such as GROMACS, LAMMPS, and Tinker-HP.

Elucidation of the noncovalent interactions driving enzyme activity guides branching enzyme engineering for α-glucan modification

Branching enzymes (BEs) confer to α-glucans, the primary energy-storage reservoir in nature, a variety of features, like slow digestion. The full catalytic cycle of BEs can be divided in six steps, namely two covalent catalytic steps involving glycosylation and transglycosylation, and four noncatalytic steps involving substrate binding and transfers (SBTs). Despite the ever-growing wealth of biochemical and structural information on BEs, clear mechanistic insights into SBTs from an industrial-performance perspective are still missing. Here, we report a Rhodothermus profundi BE (RpBE) endowed with twice as much enzymatic activity as the Rhodothermus obamensis BE currently used in industry. Furthermore, we focus on the SBTs for RpBE by means of large-scale computations supported by experiment. Engineering of the crucial positions responsible for the initial substrate-binding step improves enzymatic activity significantly, while offering a possibility to customize product types. In addition, we show that the high-efficiency substrate-transfer steps preceding glycosylation and transglycosylation are the main reason for the remarkable enzymatic activity of RpBE, suggestive of engineering directions for the BE family.

Statistical Mechanical Theories of Membrane Permeability

A popular theoretical framework to compute the permeability coefficient of a molecule is provided by the classic Smoluchowski-Kramers treatment of the steady-state diffusive flux across a free-energy barrier. Within this framework, commonly termed "inhomogeneous solubility-diffusion" (ISD), the permeability, P, is expressed in closed form in terms of the potential of mean force and position-dependent diffusivity of the molecule of interest along the membrane normal. In principle, both quantities can be calculated from all-atom MD simulations. Although several methods exist for calculating the position-dependent diffusivity, each of these is at best an estimate. In addition, the ISD model does not account for memory effects along the chosen reaction coordinate. For these reasons, it is important to seek alternative theoretical formulations to determine the permeability coefficient that are able to account for the factors ignored by the ISD approximation. Using Green-Kubo linear response theory, we establish the familiar constitutive relation between the flux density across the membrane and the difference in the concentration of a permeant molecule, j = PΔC. On this basis, we derive a time-correlation function expression for the nonequilibrium flux across a membrane that is reminiscent of the transmission coefficient in the reactive flux formalism treatment of transition rates. An analysis based on the transition path theory framework is exploited to derive alternative expressions for the permeability coefficient. The different strategies are illustrated with stochastic simulations based on the generalized Langevin equation in addition to unbiased molecular dynamics simulations of water permeation of a lipid bilayer.

Improvements in Precision of Relative Binding Free Energy Calculations Afforded by the Alchemical Enhanced Sampling (ACES) Approach

Accurate in silico predictions of how strongly small molecules bind to proteins, such as those afforded by relative binding free energy (RBFE) calculations, can greatly increase the efficiency of the hit-to-lead and lead optimization stages of the drug discovery process. The success of such calculations, however, relies heavily on their precision. Here, we show that a recently developed alchemical enhanced sampling (ACES) approach can consistently improve the precision of RBFE calculations on a large and diverse set of proteins and small molecule ligands. The addition of ACES to conventional RBFE calculations lowered the average hysteresis by over 35% (0.3-0.4 kcal/mol) and the average replicate spread by over 25% (0.2-0.3 kcal/mol) across a set of 10 protein targets and 213 small molecules while maintaining similar or improved accuracy. We show in atomic detail how ACES improved convergence of several representative RBFE calculations through enhancing the sampling of important slowly transitioning ligand degrees of freedom.

Learning Collective Variables with Synthetic Data Augmentation through Physics-Inspired Geodesic Interpolation

In molecular dynamics simulations, rare events, such as protein folding, are typically studied using enhanced sampling techniques, most of which are based on the definition of a collective variable (CV) along which acceleration occurs. Obtaining an expressive CV is crucial, but often hindered by the lack of information about the particular event, e.g., the transition from unfolded to folded conformation. We propose a simulation-free data augmentation strategy using physics-inspired metrics to generate geodesic interpolations resembling protein folding transitions, thereby improving sampling efficiency without true transition state samples. This new data can be used to improve the accuracy of classifier-based methods. Alternatively, a regression-based learning scheme for CV models can be adopted by leveraging the interpolation progress parameter.

QM/MM Free Energy Calculations of Long-Range Biological Protonation Dynamics by Adaptive and Focused Sampling

Water-mediated proton transfer reactions are central for catalytic processes in a wide range of biochemical systems, ranging from biological energy conversion to chemical transformations in the metabolism. Yet, the accurate computational treatment of such complex biochemical reactions is highly challenging and requires the application of multiscale methods, in particular hybrid quantum/classical (QM/MM) approaches combined with free energy simulations. Here, we combine the unique exploration power of new advanced sampling methods with density functional theory (DFT)-based QM/MM free energy methods for multiscale simulations of long-range protonation dynamics in biological systems. In this regard, we show that combining multiple walkers/well-tempered metadynamics with an extended system adaptive biasing force method (MWE) provides a powerful approach for exploration of water-mediated proton transfer reactions in complex biochemical systems. We compare and combine the MWE method also with QM/MM umbrella sampling and explore the sampling of the free energy landscape with both geometric (linear combination of proton transfer distances) and physical (center of excess charge) reaction coordinates and show how these affect the convergence of the potential of mean force (PMF) and the activation free energy. We find that the QM/MM-MWE method can efficiently explore both direct and water-mediated proton transfer pathways together with forward and reverse hole transfer mechanisms in the highly complex proton channel of respiratory Complex I, while the QM/MM-US approach shows a systematic convergence of selected long-range proton transfer pathways. In this regard, we show that the PMF along multiple proton transfer pathways is recovered by combining the strengths of both approaches in a QM/MM-MWE/focused US (FUS) scheme and reveals new mechanistic insight into the proton transfer principles of Complex I. Our findings provide a promising basis for the quantitative multiscale simulations of long-range proton transfer reactions in biological systems.

Lambda-ABF: Simplified, Portable, Accurate, and Cost-Effective Alchemical Free-Energy Computation

We introduce the lambda-Adaptive Biasing Force (lambda-ABF) method for the computation of alchemical free-energy differences. We propose a software implementation and showcase it on biomolecular systems. The method arises from coupling multiple-walker adaptive biasing force with λ-dynamics. The sampling of the alchemical variable is continuous and converges toward a uniform distribution, making manual optimization of the λ schedule unnecessary. Contrary to most other approaches, alchemical free-energy estimates are obtained immediately without any postprocessing. Free diffusion of λ improves orthogonal relaxation compared to fixed-λ thermodynamic integration or free-energy perturbation. Furthermore, multiple walkers provide generic orthogonal space coverage with minimal user input and negligible computational overhead. We show that our high-performance implementations coupling the Colvars library with NAMD and Tinker-HP can address real-world cases including ligand-receptor binding with both fixed-charge and polarizable models, with a demonstrably richer sampling than fixed-λ methods. The implementation is fully open-source, publicly available, and readily usable by practitioners of current alchemical methods. Thanks to the portable Colvars library, lambda-ABF presents a unified user interface regardless of the back-end (NAMD, Tinker-HP, or any software to be interfaced in the future), sparing users the effort of learning multiple interfaces. Finally, the Colvars Dashboard extension of the visual molecular dynamics (VMD) software provides an interactive monitoring and diagnostic tool for lambda-ABF simulations.

Computational Investigation of the Covalent Inhibition Mechanism of Bruton's Tyrosine Kinase by Ibrutinib

Covalent inhibitors represent a promising class of therapeutic compounds. Nonetheless, rationally designing covalent inhibitors to achieve a right balance between selectivity and reactivity remains extremely challenging. To better understand the covalent binding mechanism, a computational study is carried out using the irreversible covalent inhibitor of Bruton tyrosine kinase (BTK) ibrutinib as an example. A multi-μs classical molecular dynamics trajectory of the unlinked inhibitor is generated to explore the fluctuations of the compound associated with the kinase binding pocket. Then, the reaction pathway leading to the formation of the covalent bond with the cysteine residue at position 481 via a Michael addition is determined using the string method in collective variables on the basis of hybrid quantum mechanical-molecular mechanical (QM/MM) simulations. The reaction pathway shows a strong correlation between the covalent bond formation and the protonation/deprotonation events taking place sequentially in the covalent inhibition reaction, consistent with a 3-step reaction with transient thiolate and enolates intermediate states. Two possible atomistic mechanisms affecting deprotonation/protonation events from the thiolate to the enolate intermediate were observed: a highly correlated direct pathway involving proton transfer to the Cα of the acrylamide warhead from the cysteine involving one or a few water molecules and a more indirect pathway involving a long-lived enolate intermediate state following the escape of the proton to the bulk solution. The results are compared with experiments by simulating the long-time kinetics of the reaction using kinetic modeling.

Understanding the Reversible Binding of a Multichain Protein-Protein Complex through Free-Energy Calculations

We demonstrate that the binding affinity of a multichain protein-protein complex, insulin dimer, can be accurately predicted using a streamlined route of standard binding free-energy calculations. We find that chains A and C, which do not interact directly during binding, stabilize the insulin monomer structures and reduce the binding affinity of the two monomers, therefore enabling their reversible association. Notably, we confirm that although classical methods can estimate the binding affinity of the insulin dimer, conventional molecular dynamics, enhanced sampling algorithms, and classical geometrical routes of binding free-energy calculations may not fully capture certain aspects of the role played by the noninteracting chains in the binding dynamics. Therefore, this study not only elucidates the role of noninteracting chains in the reversible binding of the insulin dimer but also offers a methodological guide for investigating the reversible binding of multichain protein-protein complexes utilizing streamlined free-energy calculations.

Cyclization Scaffolding for Improved Vaccine Immunogen Stability: Application to Tau Protein in Alzheimer's Disease

Effective scaffolding of immunogens is crucial for generating conformationally selective antibodies through active immunization, particularly in the treatment of protein misfolding diseases such as Alzheimer's and Parkinson's disease. Previous computational work has revealed that a disorder-prone region of the tau protein, when in a stacked form, is predicted to structurally resemble a small, soluble protofibril, having conformational properties similar to those of experimental in vitro tau oligomers. Such an oligomeric structural mimic has the potential to serve as a vaccine immunogen design for Alzheimer's disease. In this study, we developed a cyclization scaffolding method in Rosetta, in which multiple cyclic peptides are stacked into a protofibril. Cyclization results in significant stabilization of protofibril-like structures by constraining the conformational space. Applying this method to the disorder-prone region of the tau fibril, we evaluated the metastability of the cyclized tau immunogen using molecular dynamics simulations, and we identified sequences of two cyclic constructs having high metastability in the protofibril. We then assessed their thermodynamic stability by computing the free energy required to separate a distal chain from the rest of the stacked structure. Our computational results, based on molecular dynamics simulations and free energy calculations, demonstrate that two cyclized constructs, cyclo-(VKSEKLDFKDRVQSKIFyN) and cyclo-(VKSEKLDFKDRVQSKIYvG) (lowercase letters indicate d-form amino acids), possess significantly increased thermodynamic stability in the protofibril over an uncyclized linear construct VKSEKLDFKDRVQSKI. The cyclization scaffolding approach proposed here holds promise as a means to effectively design immunogens for protein misfolding diseases, particularly those involving liposome-conjugated peptide constructs.

Protein ensemble modeling and analysis with MMMx

Proteins, especially of eukaryotes, often have disordered domains and may contain multiple folded domains whose relative spatial arrangement is distributed. The MMMx ensemble modeling and analysis toolbox (https://github.com/gjeschke/MMMx) can support the design of experiments to characterize the distributed structure of such proteins, starting from AlphaFold2 predictions or folded domain structures. Weak order can be analyzed with reference to a random coil model or to peptide chains that match the residue-specific Ramachandran angle distribution of the loop regions and are otherwise unrestrained. The deviation of the mean square end-to-end distance of chain sections from their average over sections of the same sequence length reveals localized compaction or expansion of the chain. The shape sampled by disordered chains is visualized by superposition in the principal axes frame of their inertia tensor. Ensembles of different sizes and with weighted conformers can be compared based on a similarity parameter that abstracts from the ensemble width.

Perspectives on Computational Enzyme Modeling: From Mechanisms to Design and Drug Development

Understanding enzyme mechanisms is essential for unraveling the complex molecular machinery of life. In this review, we survey the field of computational enzymology, highlighting key principles governing enzyme mechanisms and discussing ongoing challenges and promising advances. Over the years, computer simulations have become indispensable in the study of enzyme mechanisms, with the integration of experimental and computational exploration now established as a holistic approach to gain deep insights into enzymatic catalysis. Numerous studies have demonstrated the power of computer simulations in characterizing reaction pathways, transition states, substrate selectivity, product distribution, and dynamic conformational changes for various enzymes. Nevertheless, significant challenges remain in investigating the mechanisms of complex multistep reactions, large-scale conformational changes, and allosteric regulation. Beyond mechanistic studies, computational enzyme modeling has emerged as an essential tool for computer-aided enzyme design and the rational discovery of covalent drugs for targeted therapies. Overall, enzyme design/engineering and covalent drug development can greatly benefit from our understanding of the detailed mechanisms of enzymes, such as protein dynamics, entropy contributions, and allostery, as revealed by computational studies. Such a convergence of different research approaches is expected to continue, creating synergies in enzyme research. This review, by outlining the ever-expanding field of enzyme research, aims to provide guidance for future research directions and facilitate new developments in this important and evolving field.

Collective Variable-Based Enhanced Sampling: From Human Learning to Machine Learning

Enhanced-sampling algorithms relying on collective variables (CVs) are extensively employed to study complex (bio)chemical processes that are not amenable to brute-force molecular simulations. The selection of appropriate CVs characterizing the slow movement modes is of paramount importance for reliable and efficient enhanced-sampling simulations. In this Perspective, we first review the application and limitations of CVs obtained from chemical and geometrical intuition. We also introduce path-sampling algorithms, which can identify path-like CVs in a high-dimensional free-energy space. Machine-learning algorithms offer a viable approach to finding suitable CVs by analyzing trajectories from preliminary simulations. We discuss both the performance of machine-learning-derived CVs in enhanced-sampling simulations of experimental models and the challenges involved in applying these CVs to realistic, complex molecular assemblies. Moreover, we provide a prospective view of the potential advancements of machine-learning algorithms for the development of CVs in the field of enhanced-sampling simulations.

Enhanced Sampling of Buried Charges in Free Energy Calculations Using Replica Exchange with Charge Tempering

Buried ionizable groups in proteins often play important structural and functional roles. However, it is generally challenging to study the detailed molecular mechanisms solely based on experimental measurements. Free energy calculations using atomistic simulations, on the other hand, complement experimental studies and can provide high temporal and spatial resolution information that can lead to mechanistic insights. Nevertheless, it is also well recognized that sufficient sampling of such atomistic simulations can be challenging, considering that structural changes related to the buried charges may be very slow. In the present study, we describe a simple but effective enhanced sampling technique called replica exchange with charge tempering (REChgT) with a novel free energy method, multisite λ dynamics (MSλD), to study two systems containing buried charges, pKa prediction of a small molecule, orotate, in complex with the dihydroorotate dehydrogenase, and relative stability of a Glu-Lys pair buried in the hydrophobic core of two variants of Staphylococcal nuclease. Compared to the original MSλD simulations, the usage of REChgT dramatically increases sampling in both conformational and alchemical spaces, which directly translates into a significant reduction of wall time to converge the free energy calculations. This study highlights the importance of sufficient sampling toward developing improved free energy methods.

Classifying Protein-Protein Binding Affinity with Free-Energy Calculations and Machine Learning Approaches

Understanding the intricate phenomenon of neuronal wiring in the brain is of great interest in neuroscience. In the fruit fly, Drosophila melanogaster, the Dpr-DIP interactome has been identified to play an important role in this process. However, experimental data suggest that a merely limited subset of complexes, essentially 57 out of a total of 231, exhibit strong binding affinity. In this work, we sought to identify the residue-level molecular basis underlying the difference in binding affinity using a state-of-the-art methodology consisting of standard binding free-energy calculations with a geometrical route and machine learning (ML) techniques. We determined the binding affinity for two complexes using statistical mechanics simulations, achieving an excellent reproduction of the experimental data. Moreover, we predicted the binding free energy for two additional low-affinity complexes, devoid of experimental estimation, while simultaneously identifying key residues for the binding. Furthermore, through the use of ML algorithms, linear discriminant analysis, and random forest, we achieved remarkable accuracy, as high as 0.99, in discerning between strong (cognate) and weak (noncognate) binders. The presented ML approach encompasses easily transferable input features, enabling its broad application to any interactome while facilitating the identification of pivotal residues critical for binding interactions. The predictive power of the generated model was probed on similar protein families from 13 diverse species. Our ML model exhibited commendable performance on these additional data sets, showcasing its reliability and robustness across the species barrier.

How is Membrane Permeation of Small Ionizable Molecules Affected by Protonation Kinetics?

According to the pH-partition hypothesis, the aqueous solution adjacent to a membrane is a mixture of the ionization states of the permeating molecule at fixed Henderson-Hasselbalch concentrations, such that each state passes through the membrane in parallel with its own specific permeability. An alternative view, based on the assumption that the rate of switching ionization states is instantaneous, represents the permeation of ionizable molecules via an effective Boltzmann-weighted average potential (BWAP). Such an assumption is used in constant-pH molecular dynamics simulations. The inhomogeneous solubility-diffusion framework can be used to compute the pH-dependent membrane permeability for each of these two limiting treatments. With biased WTM-eABF molecular dynamics simulations, we computed the potential of mean force and diffusivity of each ionization state of two weakly basic small molecules: nicotine, an addictive drug, and varenicline, a therapeutic for treating nicotine addiction. At pH = 7, the BWAP effective permeability is greater than that determined by pH-partitioning by a factor of 2.5 for nicotine and 5 for varenicline. To assess the importance of ionization kinetics, we present a Smoluchowski master equation that includes explicitly the protonation and deprotonation processes coupled with the diffusive motion across the membrane. At pH = 7, the increase in permeability due to the explicit ionization kinetics is negligible for both nicotine and varenicline. This finding is reaffirmed by combined Brownian dynamics and Markov state model simulations for estimating the permeability of nicotine while allowing changes in its ionization state. We conclude that for these molecules the pH-partition hypothesis correctly captures the physics of the permeation process. The small free energy barriers for the permeation of nicotine and varenicline in their deprotonated neutral forms play a crucial role in establishing the validity of the pH-partitioning mechanism. Essentially, BWAP fails because ionization kinetics are too slow on the time scale of membrane crossing to affect the permeation of small ionizable molecules such as nicotine and varenicline. For the singly protonated state of nicotine, the computational results agree well with experimental measurements (P1 = 1.29 × 10-7 cm/s), but the agreement for neutral (P0 = 6.12 cm/s) and doubly protonated nicotine (P2 = 3.70 × 10-13 cm/s) is slightly worse, likely due to factors associated with the aqueous boundary layer (neutral form) or leaks through paracellular pathways (doubly protonated form).

Analyzing Molecular Dynamics Trajectories Thermodynamically through Artificial Intelligence

Molecular dynamics simulations produce trajectories that correspond to vast amounts of structure when exploring biochemical processes. Extracting valuable information, e.g., important intermediate states and collective variables (CVs) that describe the major movement modes, from molecular trajectories to understand the underlying mechanisms of biological processes presents a significant challenge. To achieve this goal, we introduce a deep learning approach, coined DIKI (deep identification of key intermediates), to determine low-dimensional CVs distinguishing key intermediate conformations without a-priori assumptions. DIKI dynamically plans the distribution of latent space and groups together similar conformations within the same cluster. Moreover, by incorporating two user-defined parameters, namely, coarse focus knob and fine focus knob, to help identify conformations with low free energy and differentiate the subtle distinctions among these conformations, resolution-tunable clustering was achieved. Furthermore, the integration of DIKI with a path-finding algorithm contributes to the identification of crucial intermediates along the lowest free-energy pathway. We postulate that DIKI is a robust and flexible tool that can find widespread applications in the analysis of complex biochemical processes.

Elucidating the Molecular Determinants of the Binding Modes of a Third-Generation HIV-1 Integrase Strand Transfer Inhibitor: The Importance of Side Chain and Solvent Reorganization

The first- and second-generation clinically used HIV-1 integrase (IN) strand transfer inhibitors (INSTIs) are key components of antiretroviral therapy (ART), which work by blocking the integration step in the HIV-1 replication cycle that is catalyzed by a nucleoprotein assembly called an intasome. However, resistance to even the latest clinically used INSTIs is beginning to emerge. Developmental third-generation INSTIs, based on naphthyridine scaffolds, are promising candidates to combat drug-resistant viral variants. Among these novel INSTIs, compound 4f exhibits two distinct conformations when binding with intasomes from HIV-1 and the closely related prototype foamy virus (PFV) despite the high structural similarity of their INSTI binding pockets. The molecular mechanism and the key active site residues responsible for these differing binding modes in closely related intasomes remain elusive. To unravel the molecular determinants governing the two distinct binding modes, we applied a novel molecular dynamics-based free energy method that utilizes alchemical pathways to overcome the sampling challenges associated with transitioning between the two bound conformations of ligand 4f within the crowded environments of the INSTI binding pockets in these intasomes. The calculated conformational free energies successfully recapitulate the experimentally observed binding mode preferences in the two viral intasomes. Analysis of the simulated structures suggests that the observed binding mode preferences are caused by amino acid residue differences in both the front and the central catalytic sub-pocket of the INSTI binding site in HIV-1 and PFV. Additional free energy calculations on mutants of HIV-1 and PFV revealed that while both sub-pockets contribute to binding mode selection, the central sub-pocket plays a more important role. These results highlight the importance of both side chain and solvent reorganization, as well as the conformational entropy in determining the ligand binding mode, and will help inform the development of more effective INSTIs for combatting drug-resistant viral variants.

A Rigorous Framework for Calculating Protein-Protein Binding Affinities in Membranes

Calculating the binding free energy of integral transmembrane (TM) proteins is crucial for understanding the mechanisms by which they recognize one another and reversibly associate. The glycophorin A (GpA) homodimer, composed of two α-helical segments, has long served as a model system for studying TM protein reversible association. The present work establishes a methodological framework for calculating the binding affinity of the GpA homodimer in the heterogeneous environment of a membrane. Our investigation carefully considered a variety of protocols, including the appropriate choice of the force field, rigorous standardization reflecting the experimental conditions, sampling algorithm, anisotropic environment, and collective variables, to accurately describe GpA dimerization via molecular dynamics-based approaches. Specifically, two strategies were explored: (i) an unrestrained potential mean force (PMF) calculation, which merely enhances sampling along the separation of the two binding partners without any restraint, and (ii) a so-called "geometrical route", whereby the α-helices are progressively separated with imposed restraints on their orientational, positional, and conformational degrees of freedom to accelerate convergence. Our simulations reveal that the simplified, unrestrained PMF approach is inadequate for the description of GpA dimerization. Instead, the geometrical route, tailored specifically to GpA in a membrane environment, yields excellent agreement with experimental data within a reasonable computational time. A dimerization free energy of -10.7 kcal/mol is obtained, in fairly good agreement with available experimental data. The geometrical route further helps elucidate how environmental forces drive association before helical interactions stabilize it. Our simulations also brought to light a distinct, long-lived spatial arrangement that potentially serves as an intermediate state during dimer formation. The methodological advances in the generalized geometrical route provide a powerful tool for accurate and efficient binding-affinity calculations of intricate TM protein complexes in inhomogeneous environments.

Improved Sampling of Adaptive Path Collective Variables by Stabilized Extended-System Dynamics

Because of the complicated multistep nature of many biocatalytic reactions, an a priori definition of reaction coordinates is difficult. Therefore, we apply enhanced sampling algorithms along with adaptive path collective variables (PCVs), which converge to the minimum free energy path (MFEP) during the simulation. We show how PCVs can be combined with the highly efficient well-tempered metadynamics extended-system adaptive biasing force (WTM-eABF) hybrid sampling algorithm, offering dramatically increased sampling efficiency due to its fast adaptation to path updates. For this purpose, we address discontinuities of PCVs that can arise due to path shortcutting or path updates with a novel stabilization algorithm for extended-system methods. In addition, we show how the convergence of simulations can be further accelerated by utilizing the multistate Bennett's acceptance ratio (MBAR) estimator. These methods are applied to the first step of the enzymatic reaction mechanism of pseudouridine synthases, where the ability of path WTM-eABF to efficiently explore intricate molecular transitions is demonstrated.

Binding Thermodynamics of Fourth-Generation EGFR Inhibitors Revealed by Absolute Binding Free Energy Calculations

The overexpression or mutation of the kinase domain of the epidermal growth factor receptor (EGFR) is strongly associated with non-small-cell lung cancer (NSCLC). EGFR tyrosine kinase inhibitors (TKIs) have proven to be effective in treating NSCLC patients. However, EGFR mutations can result in drug resistance. To elucidate the mechanisms underlying this resistance and inform future drug development, we examined the binding affinities of BLU-945, a recently reported fourth-generation TKI, to wild-type EGFR (EGFRWT) and its double-mutant (L858R/T790M; EGFRDM) and triple-mutant (L858R/T790M/C797S; EGFRTM) forms. We compared the binding affinities of BLU-945, BLU-945 analogues, CH7233163 (another fourth-generation TKI), and erlotinib (a first-generation TKI) using absolute binding free energy calculations. Our findings reveal that BLU-945 and CH7233163 exhibit binding affinities to both EGFRDM and EGFRTM stronger than those of erlotinib, corroborating experimental data. We identified K745 and T854 as the key residues in the binding of fourth-generation EGFR TKIs. Electrostatic forces were the predominant driving force for the binding of fourth-generation TKIs to EGFR mutants. Furthermore, we discovered that the incorporation of piperidinol and sulfone groups in BLU-945 substantially enhanced its binding capacity to EGFR mutants. Our study offers valuable theoretical insights for optimizing fourth-generation EGFR TKIs.

Standard Binding Free-Energy Calculations: How Far Are We from Automation?

Recent success stories suggest that in silico protein-ligand binding free-energy calculations are approaching chemical accuracy. However, their widespread application remains limited by the extensive human intervention required, posing challenges for the neophyte. As such, it is critical to develop automated workflows for estimating protein-ligand binding affinities with minimum personal involvement. Key human efforts include setting up and tuning enhanced-sampling or alchemical-transformation algorithms as a preamble to computational binding free-energy estimations. Additionally, preparing input files, bookkeeping, and postprocessing represent nontrivial tasks. In this Perspective, we discuss recent progress in automating standard binding free-energy calculations, featuring the development of adaptive or parameter-free algorithms, standardization of binding free-energy calculation workflows, and the implementation of user-friendly software. We also assess the current state of automated standard binding free-energy calculations and evaluate the limitations of existing methods. Last, we outline the requirements for future algorithms and workflows to facilitate automated free-energy calculations for diverse protein-ligand complexes.

Elucidating the molecular determinants for binding modes of a third-generation HIV-1 integrase strand transfer inhibitor: Importance of side chain and solvent reorganization

The first and second-generation clinically used HIV-1 integrase (IN) strand transfer inhibitors (INSTIs) are key components of antiretroviral therapy (ART), which work by blocking the integration step in the HIV-1 replication cycle that is catalyzed by a nucleoprotein assembly called an intasome. However, resistance to even the latest clinically used INSTIs is beginning to emerge. Developmental third-generation INSTIs, based on naphthyridine scaffold, are promising candidates to combat drug-resistant viral variants. Among these novel INSTIs, compound 4f exhibits two distinct conformations when binding to intasomes from HIV-1 and the closely related prototype foamy virus (PFV), despite the high structural similarity of their INSTI binding pockets. The molecular mechanism and the key active site residues responsible for these differing binding modes in closely related intasomes remain elusive. To unravel the molecular determinants governing the two distinct binding modes, we employ a novel molecular dynamics-based free energy approach that utilizes alchemical pathways to overcome the sampling challenges associated with transitioning between two ligand conformations within crowded environments along physical pathways. The calculated conformational free energies successfully recapitulate the experimentally observed binding mode preferences in the two viral intasomes. Analysis of the simulated structures suggests that the observed binding mode preferences are caused by amino acid residue differences in both the front and the central catalytic sub-pocket of the INSTI binding site in HIV-1 and PFV. Additional free energy calculations on mutants of HIV-1 and PFV revealed that while both sub-pockets contribute to the binding mode selection, the central sub-pocket plays a more important role. These results highlight the importance of both side chain and solvent reorganization, as well as the conformational entropy in determining the ligand binding mode and will help inform the development of more effective INSTIs for combatting drug-resistant viral variants.

Factorization in molecular modeling and belief propagation algorithms

Factorization reduces computational complexity, and is therefore an important tool in statistical machine learning of high dimensional systems. Conventional molecular modeling, including molecular dynamics and Monte Carlo simulations of molecular systems, is a large research field based on approximate factorization of molecular interactions. Recently, the local distribution theory was proposed to factorize joint distribution of a given molecular system into trainable local distributions. Belief propagation algorithms are a family of exact factorization algorithms for (junction) trees, and are extended to approximate loopy belief propagation algorithms for graphs with loops. Despite the fact that factorization of probability distribution is the common foundation, computational research in molecular systems and machine learning studies utilizing belief propagation algorithms have been carried out independently with respective track of algorithm development. The connection and differences among these factorization algorithms are briefly presented in this perspective, with the hope to intrigue further development of factorization algorithms for physical modeling of complex molecular systems.

Deep-Learning-Assisted Enhanced Sampling for Exploring Molecular Conformational Changes

We present a novel strategy to explore conformational changes and identify stable states of molecular objects, eliminating the need for a priori knowledge. The approach applies a deep learning method to extract information about the movement modes of the molecular object from a short, high-dimensional, and parameter-free preliminary enhanced-sampling simulation. The gathered information is described by a small set of deep-learning-based collective variables (dCVs), which steer the production-enhanced-sampling simulation. Considering the challenge of adequately exploring the configurational space using the low-dimensional, suboptimal dCVs, we incorporate a method designed for ergodic sampling, namely, Gaussian-accelerated molecular dynamics (MD), into the framework of CV-based enhanced sampling. MD simulations on both toy models and nontrivial examples demonstrate the remarkable computational efficiency of the strategy in capturing the conformational changes of molecular objects without a priori knowledge. Specifically, we achieved the blind folding of two fast folders, chignolin and villin, within a time scale of hundreds of nanoseconds and successfully reconstructed the free-energy landscapes that characterize their reversible folding. All in all, the presented strategy holds significant promise for investigating conformational changes in macromolecules, and it is anticipated to find extensive applications in the fields of chemistry and biology.

Predicting adsorption of organic compounds onto graphene and black phosphorus by molecular dynamics and machine learning

With an increase in production and application of various engineering nanomaterials (ENMs), they will inevitably be released into the environment. Adsorption of various organic chemicals onto ENMs will impact on their environmental behavior and toxicology. It is unrealistic to experimentally determine adsorption equilibrium constants (K) for the vast number of organics and ENMs due to high cost in expenditure and time. Herein, appropriate molecular dynamics (MD) methods were evaluated and selected by comparing experimental K values of seven organics adsorbed onto graphene with the MD-calculated ones. Machine learning (ML) models on K of organics adsorption onto graphene and black phosphorus nanomaterials were constructed based on a benchmark data set from the MD simulations. Lasso models based on Mordred descriptors outperformed ML models built by support vector machine, random forest, k-nearest neighbor, and gradient boosting decision tree, in terms of cross-validation coefficients (Q2 > 0.90). The Lasso models also outperformed conventional poly-parameter linear free energy relationship models for predicting logK. Compared with previous models, the Lasso models considered more compounds with different functional groups and thus have broader applicability domains. This study provides a promising way to fill the data gap in logK for chemicals adsorbed onto the ENMs.

Targeted Protein Degradation: Advances, Challenges, and Prospects for Computational Methods

The therapeutic approach of targeted protein degradation (TPD) is gaining momentum due to its potentially superior effects compared with protein inhibition. Recent advancements in the biotech and pharmaceutical sectors have led to the development of compounds that are currently in human trials, with some showing promising clinical results. However, the use of computational tools in TPD is still limited, as it has distinct characteristics compared with traditional computational drug design methods. TPD involves creating a ternary structure (protein-degrader-ligase) responsible for the biological function, such as ubiquitination and subsequent proteasomal degradation, which depends on the spatial orientation of the protein of interest (POI) relative to E2-loaded ubiquitin. Modeling this structure necessitates a unique blend of tools initially developed for small molecules (e.g., docking) and biologics (e.g., protein-protein interaction modeling). Additionally, degrader molecules, particularly heterobifunctional degraders, are generally larger than conventional small molecule drugs, leading to challenges in determining drug-like properties like solubility and permeability. Furthermore, the catalytic nature of TPD makes occupancy-based modeling insufficient. TPD consists of multiple interconnected yet distinct steps, such as POI binding, E3 ligase binding, ternary structure interactions, ubiquitination, and degradation, along with traditional small molecule properties. A comprehensive set of tools is needed to address the dynamic nature of the induced proximity ternary complex and its implications for ubiquitination. In this Perspective, we discuss the current state of computational tools for TPD. We start by describing the series of steps involved in the degradation process and the experimental methods used to characterize them. Then, we delve into a detailed analysis of the computational tools employed in TPD. We also present an integrative approach that has proven successful for degrader design and its impact on project decisions. Finally, we examine the future prospects of computational methods in TPD and the areas with the greatest potential for impact.

Entropy and Energy Profiles of Chemical Reactions

The description of chemical processes at the molecular level is often facilitated by the use of reaction coordinates or collective variables (CVs). The CV measures the progress of the reaction and allows the construction of profiles that track how specific properties evolve as the reaction progresses. Whereas CVs are routinely used, especially alongside enhanced sampling techniques, the links among reaction profiles, thermodynamic state functions, and reaction rate constants are not rigorously exploited. Here, we report a unified treatment of such reaction profiles. Tractable expressions are derived for the free-energy, internal-energy, and entropy profiles as functions of only the CV. We demonstrate the ability of this treatment to extract quantitative insight from the entropy and internal-energy profiles of various real-world physicochemical processes, including intramolecular organic reactions, ionic transport in superionic electrolytes, and molecular transport in nanoporous materials.

Recent advances on molecular dynamics-based techniques to address drug membrane permeability with atomistic detail

Several factors affect the passive membrane permeation of small molecules, including size, charge, pH, or the presence of specific chemical groups. Understanding these features is paramount to identifying or designing drug candidates with optimal ADMET properties and this can be achieved through experimental/knowledge-based methodologies or using computational approaches. Empirical methods often lack detailed information about the underlying molecular mechanism. In contrast, Molecular Dynamics-based approaches are a powerful strategy, providing an atomistic description of this process. This technique is continuously growing, featuring new related methodologies. In this work, the recent advances in this research area will be discussed.

Predictions from First-Principles of Membrane Permeability to Small Molecules: How Useful Are They in Practice?

Predicting from first-principles the rate of passive permeation of small molecules across the biological membrane represents a promising strategy for screening lead compounds upstream in the drug-discovery and development pipeline. One popular avenue for the estimation of permeation rates rests on computer simulations in conjunction with the inhomogeneous solubility-diffusion model, which requires the determination of the free-energy change and position-dependent diffusivity of the substrate along the translocation pathway through the lipid bilayer. In this Perspective, we will clarify the physical meaning of the membrane permeability inferred from such computer simulations, and how theoretical predictions actually relate to what is commonly measured experimentally. We will also examine why these calculations remain both technically challenging and overly computationally expensive, which has hitherto precluded their routine use in nonacademic settings. We finally synopsize possible research directions to meet these challenges, increase the predictive power of physics-based rates of passive permeation, and, by ricochet, improve their practical usefulness.

Discovering Reaction Pathways, Slow Variables, and Committor Probabilities with Machine Learning

A significant challenge faced by atomistic simulations is the difficulty, and often impossibility, to sample the transitions between metastable states of the free-energy landscape associated with slow molecular processes. Importance-sampling schemes represent an appealing option to accelerate the underlying dynamics by smoothing out the relevant free-energy barriers, but require the definition of suitable reaction-coordinate (RC) models expressed in terms of compact low-dimensional sets of collective variables (CVs). While most computational studies of slow molecular processes have traditionally relied on educated guesses based on human intuition to reduce the dimensionality of the problem at hand, a variety of machine-learning (ML) algorithms have recently emerged as powerful alternatives to discover meaningful CVs capable of capturing the dynamics of the slowest degrees of freedom. Considering a simple paradigmatic situation in which the long-time dynamics is dominated by the transition between two known metastable states, we compare two variational data-driven ML methods based on Siamese neural networks aimed at discovering a meaningful RC model─the slowest decorrelating CV of the molecular process, and the committor probability to first reach one of the two metastable states. One method is the state-free reversible variational approach for Markov processes networks (VAMPnets), or SRVs─the other, inspired by the transition path theory framework, is the variational committor-based neural networks, or VCNs. The relationship and the ability of these methodologies to discover the relevant descriptors of the slow molecular process of interest are illustrated with a series of simple model systems. We also show that both strategies are amenable to importance-sampling schemes through an appropriate reweighting algorithm that approximates the kinetic properties of the transition.

A Consistent Thermodynamic Characterization of the Adsorption Process through the Calculation of Free Energy Contributions

We apply in this study different methodologies based on thermodynamic integration (TI), free energy perturbation (FEP), and potential of mean force (PMF) to address the challenging issue of the calculation of the free energy of adsorption. A model system composed of a solid substrate, an adsorbate, and solvent particles is specifically designed to reduce the dependence of our free energy results on the sampling of the phase space and the choice of the pathway. The reliability and efficiency of these alchemical free energy simulations are established through the closure of a thermodynamic cycle describing the adsorption process in solution and in a vacuum. We complete this study by the calculation of free energy contributions related to phenomena of desorption of solvent molecules and desolvation of the adsorbate upon adsorption. This calculation relies on the work of adhesion, the interfacial tension of the liquid-vapor of the solvent, and the free energy of solvation of the substrate. The different ways of calculating the free energy of adsorption are in excellent agreement and could complete experiments in the field of adsorption by giving quantitative data on the different energy contributions involved in the process.

Improving Speed and Affordability without Compromising Accuracy: Standard Binding Free-Energy Calculations Using an Enhanced Sampling Algorithm, Multiple-Time Stepping, and Hydrogen Mass Repartitioning

Accurate evaluation of protein-ligand binding free energies in silico is of paramount importance for understanding the mechanisms of biological regulation and providing a theoretical basis for drug design and discovery. Based on a series of atomistic molecular dynamics simulations in an explicit solvent, using well-tempered metadynamics extended adaptive biasing force (WTM-eABF) as an enhanced sampling algorithm, the so-called "geometrical route" offers a rigorous theoretical framework for binding affinity calculations that match experimental values. However, although robust, this strategy remains expensive, requiring substantial computational time to achieve convergence of the simulations. Improving the efficiency of the geometrical route, while preserving its reliability through improved ergodic sampling, is, therefore, highly desirable. In this contribution, having identified the computational bottleneck of the geometrical route, to accelerate the calculations we combine (i) a longer time step for the integration of the equations of motion with hydrogen-mass repartitioning (HMR), and (ii) multiple time-stepping (MTS) for collective-variable and biasing-force evaluation. Altogether, we performed 50 independent WTM-eABF simulations in triplicate for the "physical" separation of the Abl kinase-SH3 domain:p41 complex, following different HMR and MTS schemes, while tuning, in distinct protocols, the parameters of the enhanced-sampling algorithm. To demonstrate the consistency and reliability of the results obtained with the best-performing setups, we carried out quintuple simulations. Furthermore, we demonstrated the transferability of our method to other complexes by triplicating a 200 ns separation simulation of nine chosen protocols for the MDM2-p53:NVP-CGM097 complex. [Holzer et al. J. Med. Chem. 2015, 58, 6348-6358.] Our results, based on an aggregate simulation time of 14.4 μs, allowed an optimal set of parameters to be identified, able to accelerate convergence by a factor of three without any noticeable loss of accuracy.

Sigmoid Accelerated Molecular Dynamics: An Efficient Enhanced Sampling Method for Biosystems

Gaussian accelerated molecular dynamics (GaMD) is recognized as a popular enhanced sampling method for tackling long-standing challenges in biomolecular simulations. Inspired by GaMD, Sigmoid accelerated molecular dynamics (SaMD) is proposed in this work by adding a Sigmoid boost potential to improve the balance between the highest acceleration and accurate reweighting. Compared with GaMD, SaMD extends the accessible time scale and improves the computational efficiency as tested in three tasks. In the alanine dipeptide task, SaMD can produce the free energy landscape with better accuracy and efficiency. In the chignolin folding task, the estimated Gibbs free energy difference can converge to the experimental value ∼30% faster. In the protein-ligand binding task, the bound conformations are closer to the crystal structure with a minimal ligand root-mean-square deviation of 1.7 Å. The binding of the ligand XK263 to the HIV protease is reproduced by SaMD in ∼60% less simulation time.

Challenges in Protein QM/MM Simulations with Intra-Backbone Link Atoms

Hybrid quantum mechanical/molecular mechanical (QM/MM) simulations fuel discoveries in many fields of science including computational biochemistry and enzymology. Development of more convenient tools leads to an increase in the number of works in which mechanical insights into enzymes' mode of operation are obtained. Most commonly, these tools feature hydrogen-capping (link atom) approach to provide coupling between QM and MM subsystems across a covalent bond. Extensive studies were conducted to provide a solid foundation for the correctness of such an approach when a bond to a nonpolar MM atom is considered. However, not every task may be accomplished this way. Certain scenarios of using QM/MM in computational enzymology encourage or even necessitate the incorporation of backbone atoms into the QM region. Two out of three backbone atoms are polar, and in QM/MM with electrostatic embedding, a neighboring link atom will be hyperpolarized. Several schemes to mitigate this effect were previously proposed alongside a rigorous assessment of quantitative effects on model systems. However, it was not clear whether they may translate into qualitatively different results and how link atom hyperpolarization may manifest itself in a real-life enzymological scenario. Here, we show that the consequences of such an artifact may be severe and may completely overturn the conclusions drawn from the simulations. Our case advocates for the use of charge redistribution schemes whenever intra-backbone QM/MM boundaries are considered. Moreover, we addressed how different boundary types and charge redistribution schemes influence backbone dynamics. We showed that the results are heavily dependent on which boundary MM terms are retained, with charge alteration being of secondary importance. In the worst case, only three intra-backbone boundaries may be used with relative confidence in the adequacy of resulting simulations, irrespective of the hyperpolarization mitigation scheme. Thus, advances in the field are certainly needed to fuel new discoveries. As of now, we believe that issues raised in this work might encourage authors in the field to report what boundaries, boundary MM terms, and charge redistribution schemes they are using, so their results may be correctly interpreted.

Cooperation of structural motifs controls drug selectivity in cyclin-dependent kinases: an advanced theoretical analysis

Understanding drug selectivity mechanism is a long-standing issue for helping design drugs with high specificity. Designing drugs targeting cyclin-dependent kinases (CDKs) with high selectivity is challenging because of their highly conserved binding pockets. To reveal the underlying general selectivity mechanism, we carried out comprehensive analyses from both the thermodynamics and kinetics points of view on a representative CDK12 inhibitor. To fully capture the binding features of the drug-target recognition process, we proposed to use kinetic residue energy analysis (KREA) in conjunction with the community network analysis (CNA) to reveal the underlying cooperation effect between individual residues/protein motifs to the binding/dissociating process of the ligand. The general mechanism of drug selectivity in CDKs can be summarized as that the difference of structural cooperation between the ligand and the protein motifs leads to the difference of the energetic contribution of the key residues to the ligand. The proposed mechanisms may be prevalent in drug selectivity issues, and the insights may help design new strategies to overcome/attenuate the drug selectivity associated problems.

Perturbation of Lipid Bilayers by Biomimetic Photoswitches Based on Cyclocurcumin

The use of photoswitches which may be activated by suitable electromagnetic radiation is an attractive alternative to conventional photodynamic therapy. Here, we report all-atom molecular dynamics simulation of a biomimetic photoswitch derived from cyclocurcumin and experiencing E/Z photoisomerization. In particular, we show that the two isomers interact persistently with a lipid bilayer modeling a cellular membrane. Furthermore, the interaction with the membrane is strongly dependent on the concentration, and a transition between ordered and disordered arrangements of the photoswitches is observed. We also confirm that the structural parameters of the bilayer are differently affected by the two isomers and hence can be modulated through photoswitching, offering interesting perspectives for future applications.

Chasing collective variables using temporal data-driven strategies

The convergence of free-energy calculations based on importance sampling depends heavily on the choice of collective variables (CVs), which in principle, should include the slow degrees of freedom of the biological processes to be investigated. Autoencoders (AEs), as emerging data-driven dimension reduction tools, have been utilised for discovering CVs. AEs, however, are often treated as black boxes, and what AEs actually encode during training, and whether the latent variables from encoders are suitable as CVs for further free-energy calculations remains unknown. In this contribution, we review AEs and their time-series-based variants, including time-lagged AEs (TAEs) and modified TAEs, as well as the closely related model variational approach for Markov processes networks (VAMPnets). We then show through numerical examples that AEs learn the high-variance modes instead of the slow modes. In stark contrast, time series-based models are able to capture the slow modes. Moreover, both modified TAEs with extensions from slow feature analysis and the state-free reversible VAMPnets (SRVs) can yield orthogonal multidimensional CVs. As an illustration, we employ SRVs to discover the CVs of the isomerizations of N-acetyl-N'-methylalanylamide and trialanine by iterative learning with trajectories from biased simulations. Last, through numerical experiments with anisotropic diffusion, we investigate the potential relationship of time-series-based models and committor probabilities.

Docking and molecular dynamics simulation for therapeutic repurposing in small cell lung cancer (SCLC) patients infected with COVID-19

Cancer care has become a challenge with the current COVID-19 pandemic scenario. Specially, cancers like small cell lung cancers (SCLC) are difficult to treat even in the normal situation due to their rapid growth and early metastasis. For such patients, treatment can't be compromised and care must be taken to ensure their minimum exposure to the ongoing spread of COVID-19 infection. For this reason, in-house treatments are being suggested for these patients. Another issue is that symptoms of SCLC match well with that of COVID-19 infection. Hence, the detection of COVID-19 may also get delayed leading to unnecessary complications. Thus, we have tried to investigate if the therapeutics that is currently used in lung cancer treatment can also act against SARS-CoV-2. If it is so, the same treatment protocols can be continued even if the SCLC patient had contracted COVID-19 without compromising the cancer care. For this, RNA dependent RNA polymerase (RdRP) from SARS-CoV-2 has been selected as drug target. Both docking and molecular dynamicssimulation analysis have indicated that Paclitaxel and Dacomitinib may be explored as multi-target drugs for both SCLC and COVID-19.Communicated by Ramaswamy H. Sarma.

Standard Binding Free Energy and Membrane Desorption Mechanism for a Phospholipase C

Peripheral membrane proteins (PMPs) bind temporarily to cellular membranes and play important roles in signaling, lipid metabolism, and membrane trafficking. Obtaining accurate membrane-PMP affinities using experimental techniques is more challenging than for protein-ligand affinities in an aqueous solution. At the theoretical level, calculation of the standard protein-membrane binding free energy using molecular dynamics simulations remains a daunting challenge owing to the size of the biological objects at play, the slow lipid diffusion, and the large variation in configurational entropy that accompanies the binding process. To overcome these challenges, we used a computational framework relying on a series of potential-of-mean-force (PMF) calculations including a set of geometrical restraints on collective variables. This methodology allowed us to determine the standard binding free energy of a PMP to a phospholipid bilayer using an all-atom force field. Bacillus thuringiensis phosphatidylinositol-specific phospholipase C (BtPI-PLC) was chosen due to its importance as a virulence factor and owing to the host of experimental affinity data available. We computed a standard binding free energy of -8.2 ± 1.4 kcal/mol in reasonable agreement with the reported experimental values (-6.6 ± 0.2 kcal/mol). In light of the 2.3-μs separation PMF calculation, we investigated the mechanism whereby BtPI-PLC disengages from interactions with the lipid bilayer during separation. We describe how a short amphipathic helix engages in transitory interactions to ease the passage of its hydrophobes through the interfacial region upon desorption from the bilayer.

Biomimetic Photo-Switches Softening Model Lipid Membranes

We report the synthesis and characterization of a novel photo-switch based on biomimetic cyclocurcumin analogous and interacting with the lipid bilayer, which can be used in the framework of oxygen-independent light-induced therapy. More specifically, by using molecular dynamics simulations and free energy techniques, we show that the inclusion of hydrophobic substituents is needed to allow insertion in the lipid membrane. After having confirmed experimentally that the substituents do not preclude the efficient photoisomerization, we show through UV-vis and dynamic light scattering measurements together with compression isotherms that the chromophore is internalized in both lipid vesicles and monomolecular film, respectively, inducing their fluidification. The irradiation of the chromophore-loaded lipid aggregates modifies their properties due to the different organization of the two diastereoisomers, E and Z. In particular, a competition between a fast structural reorganization and a slower expulsion of the chromophore after isomerization can be observed in the kinetic profiles recorded during E to Z photoisomerization. This report paves the way for future investigations in the optimization of biomimetic photoswitches potentially useful in modern light-induced therapeutic strategies.

Enhancing sampling with free-energy calculations

In recent years, considerable progress has been made to enhance sampling and help address biological questions, including, but not limited to conformational transitions in biomolecules and protein-ligand reversible binding, hitherto intractable by brute-force computer simulations. Many of these advances result from the development of a palette of methods aimed at exploring rare events through reliable free-energy calculations. The advent of new, often conceptually related methods has also rendered difficult the choice of the best suited option for a given problem. Here, we focus on geometrical transformations and algorithms designed to enhance sampling along adequately chosen progress variables, tracing their theoretical foundations, and showing how they are connected and can be blended together for improved performance.

Nature-Derived Compounds as Potential Bioactive Leads against CDK9-Induced Cancer: Computational and Network Pharmacology Approaches

Given the importance of cyclin-dependent kinases (CDKs) in the maintenance of cell development, gene transcription, and other essential biological operations, CDK blockers have been generated to manage a variety of disorders resulting from CDK irregularities. Furthermore, CDK9 has a crucial role in transcription by regulating short-lived anti-apoptotic genes necessary for cancer cell persistence. Addressing CDK9 with blockers has consequently emerged as a promising treatment for cancer. This study scrutinizes the effectiveness of nature-derived compounds (geniposidic acid, quercetin, geniposide, curcumin, and withanolide C) against CDK9 through computational approaches. A molecular docking study was performed after preparing the protein and the ligands. The selected blockers of the CDK9 exerted reliable binding affinities (−8.114 kcal/mol to −13.908 kcal/mol) against the selected protein, resulting in promising candidates compared to the co-crystallized ligand (LCI). The binding affinity of geniposidic acid (−13.908 kcal/mol) to CDK9 is higher than quercetin (−10.775 kcal/mol), geniposide (−9.969 kcal/mol), curcumin (−9.898 kcal/mol), withanolide C (−8.114 kcal/mol), and the co-crystallized ligand LCI (−11.425 kcal/mol). Therefore, geniposidic acid is a promising inhibitor of CDK9. Moreover, the molecular dynamics studies assessed the structure–function relationships and protein–ligand interactions. The network pharmacology study for the selected ligands demonstrated the auspicious compound–target–pathway signaling pathways vital in developing tumor, tumor cell growth, differentiation, and promoting tumor cell progression. Moreover, this study concluded by analyzing the computational approaches the natural-derived compounds that have potential interacting activities against CDK9 and, therefore, can be considered promising candidates for CKD9-induced cancer. To substantiate this study’s outcomes, in vivo research is recommended.