Gaussian Accelerated Molecular Dynamics in NAMD

GENESIS 2.1: High-Performance Molecular Dynamics Software for Enhanced Sampling and Free-Energy Calculations for Atomistic, Coarse-Grained, and Quantum Mechanics/Molecular Mechanics Models

GENeralized-Ensemble SImulation System (GENESIS) is a molecular dynamics (MD) software developed to simulate the conformational dynamics of a single biomolecule, as well as molecular interactions in large biomolecular assemblies and between multiple biomolecules in cellular environments. To achieve the latter purpose, the earlier versions of GENESIS emphasized high performance in atomistic MD simulations on massively parallel supercomputers, with or without graphics processing units (GPUs). Here, we implemented multiscale MD simulations that include atomistic, coarse-grained, and hybrid quantum mechanics/molecular mechanics (QM/MM) calculations. They demonstrate high performance and are integrated with enhanced conformational sampling algorithms and free-energy calculations without using external programs except for the QM programs. In this article, we review new functions, molecular models, and other essential features in GENESIS version 2.1 and discuss ongoing developments for future releases.

Learning QM/MM potential using equivariant multiscale model

The machine learning (ML) method emerges as an efficient and precise surrogate model for high-level electronic structure theory. Its application has been limited to closed chemical systems without considering external potentials from the surrounding environment. To address this limitation and incorporate the influence of external potentials, polarization effects, and long-range interactions between a chemical system and its environment, the first two terms of the Taylor expansion of an electrostatic operator have been used as extra input to the existing ML model to represent the electrostatic environments. However, high-order electrostatic interaction is often essential to account for external potentials from the environment. The existing models based only on invariant features cannot capture significant distribution patterns of the external potentials. Here, we propose a novel ML model that includes high-order terms of the Taylor expansion of an electrostatic operator and uses an equivariant model, which can generate a high-order tensor covariant with rotations as a base model. Therefore, we can use the multipole-expansion equation to derive a useful representation by accounting for polarization and intermolecular interaction. Moreover, to deal with long-range interactions, we follow the same strategy adopted to derive long-range interactions between a target system and its environment media. Our model achieves higher prediction accuracy and transferability among various environment media with these modifications.

Unveiling Conformational States of CDK6 Caused by Binding of Vcyclin Protein and Inhibitor by Combining Gaussian Accelerated Molecular Dynamics and Deep Learning

CDK6 plays a key role in the regulation of the cell cycle and is considered a crucial target for cancer therapy. In this work, conformational transitions of CDK6 were identified by using Gaussian accelerated molecular dynamics (GaMD), deep learning (DL), and free energy landscapes (FELs). DL finds that the binding pocket as well as the T-loop binding to the Vcyclin protein are involved in obvious differences of conformation contacts. This result suggests that the binding pocket of inhibitors (LQQ and AP9) and the binding interface of CDK6 to the Vcyclin protein play a key role in the function of CDK6. The analyses of FELs reveal that the binding pocket and the T-loop of CDK6 have disordered states. The results from principal component analysis (PCA) indicate that the binding of the Vcyclin protein affects the fluctuation behavior of the T-loop in CDK6. Our QM/MM-GBSA calculations suggest that the binding ability of LQQ to CDK6 is stronger than AP9 with or without the binding of the Vcyclin protein. Interaction networks of inhibitors with CDK6 were analyzed and the results reveal that LQQ contributes more hydrogen binding interactions (HBIs) and hot interaction spots with CDK6. In addition, the binding pocket endures flexibility changes from opening to closing states and the Vcyclin protein plays an important role in the stabilizing conformation of the T-loop. We anticipate that this work could provide useful information for further understanding the function of CDK6 and developing new promising inhibitors targeting CDK6.

Recent advances in computational and experimental protein-ligand affinity determination techniques

INTRODUCTION

Modern drug discovery revolves around designing ligands that target the chosen biomolecule, typically proteins. For this, the evaluation of affinities of putative ligands is crucial. This has given rise to a multitude of dedicated computational and experimental methods that are constantly being developed and improved.

AREAS COVERED

In this review, the authors reassess both the industry mainstays and the newest trends among the methods for protein - small-molecule affinity determination. They discuss both computational affinity predictions and experimental techniques, describing their basic principles, main limitations, and advantages. Together, this serves as initial guide to the currently most popular and cutting-edge ligand-binding assays employed in rational drug design.

EXPERT OPINION

The affinity determination methods continue to develop toward miniaturization, high-throughput, and in-cell application. Moreover, the availability of data analysis tools has been constantly increasing. Nevertheless, cross-verification of data using at least two different techniques and careful result interpretation remain of utmost importance.

Base pair dynamics, electrostatics, and thermodynamics at the LTR-III quadruplex:duplex junction

G-quadruplexes (GQs) play key regulatory roles within the human genome and have also been identified to play similar roles in other eukaryotes, bacteria, archaea, and viruses. Human immunodeficiency virus 1, the etiological agent of acquired immunodeficiency syndrome, can form two GQs in its long terminal repeat (LTR) promoter region, each of which act to regulate viral gene expression in opposing manners. The major LTR GQ, called LTR-III, is a distinct hybrid GQ containing a 12-nucleotide duplex loop attached to the quadruplex motif. The resulting quadruplex:duplex junction (QDJ) has been hypothesized to serve as a selective drug targeting site. To better understand the dynamics of this QDJ, we performed conventional and enhanced-sampling molecular dynamics simulations using the Drude-2017 force field. We observed unbiased and reversible formation of additional base pairs in the QDJ, between Ade4:Thy14 and Gua3:Thy14. Both base pairs were electrostatically favored, but geometric constraints within the junction may drive the formation of, and preference for, the Ade4:Thy14 base pair. Finally, we demonstrated that the base pairs are separated only by small energy barriers that may enable transitions between both base-paired states. Together, these simulations provide new insights into the dynamics, electrostatics, and thermodynamics of the LTR-III QDJ.

Gaussian accelerated molecular dynamics simulations facilitate prediction of the permeability of cyclic peptides

Despite their widespread use as therapeutics, clinical development of small molecule drugs remains challenging. Among the many parameters that undergo optimization during the drug development process, increasing passive cell permeability (i.e., log(P)) can have some of the largest impact on potency. Cyclic peptides (CPs) have emerged as a viable alternative to small molecules, as they retain many of the advantages of small molecules (oral availability, target specificity) while being highly effective at traversing the plasma membrane. However, the relationship between the dominant conformations that typify CPs in an aqueous versus a membrane environment and cell permeability remain poorly characterized. In this study, we have used Gaussian accelerated molecular dynamics (GaMD) simulations to characterize the effect of solvent on the free energy landscape of lariat peptides, a subset of CPs that have recently shown potential for drug development (Kelly et al., JACS 2021). Differences in the free energy of lariat peptides as a function of solvent can be used to predict permeability of these molecules, and our results show that permeability is most greatly influenced by N-methylation and exposure to solvent. Our approach lays the groundwork for using GaMD as a way to virtually screen large libraries of CPs and drive forward development of CP-based therapeutics.

Activation of Polycystin-1 Signaling by Binding of Stalk-derived Peptide Agonists

Polycystin-1 (PC1) is the membrane protein product of the PKD1 gene whose mutation is responsible for 85% of the cases of autosomal dominant polycystic kidney disease (ADPKD). ADPKD is primarily characterized by the formation of renal cysts and potential kidney failure. PC1 is an atypical G protein-coupled receptor (GPCR) consisting of 11 transmembrane helices and an autocatalytic GAIN domain that cleaves PC1 into extracellular N-terminal (NTF) and membrane-embedded C-terminal (CTF) fragments. Recently, signaling activation of the PC1 CTF was shown to be regulated by a stalk tethered agonist (TA), a distinct mechanism observed in the adhesion GPCR family. A novel allosteric activation pathway was elucidated for the PC1 CTF through a combination of Gaussian accelerated molecular dynamics (GaMD), mutagenesis and cellular signaling experiments. Here, we show that synthetic, soluble peptides with 7 to 21 residues derived from the stalk TA, in particular, peptides including the first 9 residues (p9), 17 residues (p17) and 21 residues (p21) exhibited the ability to re-activate signaling by a stalkless PC1 CTF mutant in cellular assays. To reveal molecular mechanisms of stalk peptide-mediated signaling activation, we have applied a novel Peptide GaMD (Pep-GaMD) algorithm to elucidate binding conformations of selected stalk peptide agonists p9, p17 and p21 to the stalkless PC1 CTF. The simulations revealed multiple specific binding regions of the stalk peptide agonists to the PC1 protein including an "intermediate" bound yet inactive state. Our Pep-GaMD simulation findings were consistent with the cellular assay experimental data. Binding of peptide agonists to the TOP domain of PC1 induced close TOP-putative pore loop interactions, a characteristic feature of the PC1 CTF signaling activation mechanism. Using sequence covariation analysis of PC1 homologs, we further showed that the peptide binding regions were consistent with covarying residue pairs identified between the TOP domain and the stalk TA. Therefore, structural dynamic insights into the mechanisms of PC1 activation by stalk-derived peptide agonists have enabled an in-depth understanding of PC1 signaling. They will form a foundation for development of PC1 as a therapeutic target for the treatment of ADPKD.

Mechanistic Insights into Peptide Binding and Deactivation of an Adhesion G Protein-Coupled Receptor

Adhesion G protein-coupled receptors (ADGRGs) play critical roles in the reproductive, neurological, cardiovascular, and endocrine systems. In particular, ADGRG2 plays a significant role in Ewing sarcoma cell proliferation, parathyroid cell function, and male fertility. In 2022, a cryo-EM structure was reported for the active ADGRG2 bound by an optimized peptide agonist IP15 and the Gs protein. The IP15 peptide agonist was also modified to antagonists 4PH-E and 4PH-D with mutations of the 4PH residue to Glu and Asp, respectively. However, experimental structures of inactive antagonist-bound ADGRs remain to be resolved, and the activation mechanism of ADGRs such as ADGRG2 is poorly understood. Here, we applied Gaussian accelerated molecular dynamics (GaMD) simulations to probe conformational dynamics of the agonist- and antagonist-bound ADGRG2. By performing GaMD simulations, we were able to identify important low-energy conformations of ADGRG2 in the active, intermediate, and inactive states, as well as explore the binding conformations of each peptide. Moreover, our simulations revealed critical peptide-receptor residue interactions during the deactivation of ADGRG2. In conclusion, through GaMD simulations, we uncovered mechanistic insights into peptide (agonist and antagonist) binding and deactivation of the ADGRG2. These findings will potentially facilitate rational design of new peptide modulators of ADGRG2 and other ADGRs.

Development of N-Terminally Modified Variants of the CXCR4-Antagonistic Peptide EPI-X4 for Enhanced Plasma Stability

EPI-X4, a natural peptide CXCR4 antagonist, shows potential for treating inflammation and cancer, but its short plasma stability limits its clinical application. We aimed to improve the plasma stability of EPI-X4 analogues without compromising CXCR4 antagonism. Our findings revealed that only the peptide N-terminus is prone to degradation. Consequently, incorporating d-amino acids or acetyl groups in this region enhanced peptide stability in plasma. Notably, EPI-X4 leads 5, 27, and 28 not only retained their CXCR4 binding and antagonism but also remained stable in plasma for over 8 h. Molecular dynamic simulations showed that these modified analogues bind similarly to CXCR4 as the original peptide. To further increase their systemic half-lives, we conjugated these stabilized analogues with large polymers and albumin binders. These advances highlight the potential of the optimized EPI-X4 analogues as promising CXCR4-targeted therapeutics and set the stage for more detailed preclinical assessments.

Bacterial Proteases as Potentially Exploitable Modulators of SARS-CoV-2 Infection: Logic from the Literature, Informatics, and Inspiration from the Dog

(1) Background: The COVID-19 pandemic left many intriguing mysteries. Retrospective vulnerability trends tie as strongly to odd demographics as to exposure profiles, genetics, health, or prior medical history. This article documents the importance of nasal microbiome profiles in distinguishing infection rate trends among differentially affected subgroups. (2) Hypothesis: From a detailed literature survey, microbiome profiling experiments, bioinformatics, and molecular simulations, we propose that specific commensal bacterial species in the Pseudomonadales genus confer protection against SARS-CoV-2 infections by expressing proteases that may interfere with the proteolytic priming of the Spike protein. (3) Evidence: Various reports have found elevated Moraxella fractions in the nasal microbiomes of subpopulations with higher resistance to COVID-19 (e.g., adolescents, COVID-19-resistant children, people with strong dietary diversity, and omnivorous canines) and less abundant ones in vulnerable subsets (the elderly, people with narrower diets, carnivorous cats and foxes), along with bioinformatic evidence that Moraxella bacteria express proteases with notable homology to human TMPRSS2. Simulations suggest that these proteases may proteolyze the SARS-CoV-2 spike protein in a manner that interferes with TMPRSS2 priming.

Development and Evaluation of a Self-Nanoemulsifying Drug Delivery System for Sinapic Acid with Improved Antiviral Efficacy against SARS-CoV-2

This study aimed to develop a self-nanoemulsifying drug delivery system (SNE) for sinapic acid (SA) to improve its solubility and antiviral activity. Optimal components for the SA-SNE formulation were selected, including Labrafil as the oil, Cremophor EL as the surfactant, and Transcutol as the co-surfactant. The formulation was optimized using surface response design, and the optimized SA-SNE formulation exhibited a small globule size of 83.6 nm, high solubility up to 127.1 ± 3.3, and a 100% transmittance. In vitro release studies demonstrated rapid and high SA release from the formulation. Pharmacokinetic analysis showed improved bioavailability by 2.43 times, and the optimized SA-SNE formulation exhibited potent antiviral activity against SARS-CoV-2. The developed SA-SNE formulation can enhance SA's therapeutic efficacy by improving its solubility, bioavailability, and antiviral activity. Further in silico, modeling, and Gaussian accelerated molecular dynamics (GaMD)-based studies revealed that SA could interact with and inhibit the viral main protease (Mpro). This research contributes to developing effective drug delivery systems for poorly soluble drugs like SA, opening new possibilities for their application via nebulization in SARS-CoV-2 therapy.

A drug design strategy based on molecular docking and molecular dynamics simulations applied to development of inhibitor against triple-negative breast cancer by Scutellarein derivatives

Triple-negative breast cancer (TNBC), accounting for 10-15% of all breast malignancies, is more prevalent in women under 40, particularly in those of African descent or carrying the BRCA1 mutation. TNBC is characterized by the absence of estrogen and progesterone receptors (ER, PR) and low or elevated HER2 expression. It represents a particularly aggressive form of breast cancer with limited therapeutic options and a poorer prognosis. In our study, we utilized the protein of TNBC collected from the Protein Data Bank (PDB) with the most stable configuration. We selected Scutellarein, a bioactive molecule renowned for its anti-cancer properties, and used its derivatives to design potential anti-cancer drugs employing computational tools. We applied and modified structural activity relationship methods to these derivatives and evaluated the probability of active (Pa) and inactive (Pi) outcomes using pass prediction scores. Furthermore, we employed in-silico approaches such as the assessment of absorption, distribution, metabolism, excretion, and toxicity (ADMET) parameters, and quantum calculations through density functional theory (DFT). Within the DFT calculations, we analyzed Frontier Molecular Orbitals, specifically the Highest Occupied Molecular Orbital (HOMO) and Lowest Unoccupied Molecular Orbital (LUMO). We then conducted molecular docking and dynamics against TNBC to ascertain binding affinity and stability. Our findings indicated that Scutellarein derivatives, specifically DM03 with a binding energy of -10.7 kcal/mol and DM04 with -11.0 kcal/mol, exhibited the maximum binding tendency against Human CK2 alpha kinase (PDB ID 7L1X). Molecular dynamic simulations were performed for 100 ns, and stability was assessed using root-mean-square deviation (RMSD) and root-mean-square fluctuation (RMSF) parameters, suggesting significant stability for our chosen compounds. Furthermore, these molecules met the pharmacokinetics requirements for potential therapeutic candidates, displaying non-carcinogenicity, minimal aquatic and non-aquatic toxicity, and greater aqueous solubility. Collectively, our computational data suggest that Scutellarein derivatives may serve as potential therapeutic agents for TNBC. However, further experimental investigations are needed to validate these findings.

Flexible Gaussian Accelerated Molecular Dynamics to Enhance Biological Sampling

Molecular dynamics simulations often struggle to obtain sufficient sampling to study complex molecular events due to high energy barriers separating the minima of interest. Multiple enhanced sampling techniques have been developed and improved over the years to tackle this issue. Gaussian accelerated molecular dynamics (GaMD) is a recently developed enhanced sampling technique that works by adding a biasing potential, lifting the energy landscape up, and decreasing the height of its barriers. GaMD allows one to increase the sampling of events of interest without the need of a priori knowledge of the system or the relevant coordinates. All required acceleration parameters can be obtained from a previous search run. Upon its development, several improvements for the methodology have been proposed, among them selective GaMD in which the boosting potential is selectively applied to the region of interest. There are currently four selective GaMD methods that have shown promising results. However, all of these methods are constrained on the number, location, and scenarios in which this selective boosting potential can be applied to ligands, peptides, or protein-protein interactions. In this work, we showcase a GROMOS implementation of the GaMD methodology with a fully flexible selective GaMD approach that allows the user to define, in a straightforward way, multiple boosting potentials for as many regions as desired. We show and analyze the advantages of this flexible selective approach on two previously used test systems, the alanine dipeptide and the chignolin peptide, and extend these examples to study its applicability and potential to study conformational changes of glycans and glycosylated proteins.

Temperature Accelerated Sliced Sampling to Probe Ligand Dissociation from Protein

Modeling ligand unbinding in proteins to estimate the free energy of binding and probing the mechanism presents several challenges. They primarily pertain to the entropic bottlenecks resulting from protein and solvent conformations. While exploring the unbinding processes using enhanced sampling techniques, very long simulations are required to sample all of the conformational states as the system gets trapped in local free energy minima along transverse coordinates. Here, we demonstrate that temperature accelerated sliced sampling (TASS) is an ideal approach to overcome some of the difficulties faced by conventional sampling methods in studying ligand unbinding. Using TASS, we study the unbinding of avibactam inhibitor molecules from the Class C β-lactamase (CBL) active site. Extracting CBL-avibactam unbinding free energetics, unbinding pathways, and identifying critical interactions from the TASS simulations are demonstrated.

Spermine and spermidine bind CXCR4 and inhibit CXCR4- but not CCR5-tropic HIV-1 infection

Semen is an important vector for sexual HIV-1 transmission. Although CXCR4-tropic (X4) HIV-1 may be present in semen, almost exclusively CCR5-tropic (R5) HIV-1 causes systemic infection after sexual intercourse. To identify factors that may limit sexual X4-HIV-1 transmission, we generated a seminal fluid-derived compound library and screened it for antiviral agents. We identified four adjacent fractions that blocked X4-HIV-1 but not R5-HIV-1 and found that they all contained spermine and spermidine, abundant polyamines in semen. We showed that spermine, which is present in semen at concentrations up to 14 mM, binds CXCR4 and selectively inhibits cell-free and cell-associated X4-HIV-1 infection of cell lines and primary target cells at micromolar concentrations. Our findings suggest that seminal spermine restricts sexual X4-HIV-1 transmission.

Multivalent Molecular Tweezers Disrupt the Essential NDC80 Interaction with Microtubules

Binding of microtubule filaments by the conserved Ndc80 protein is required for kinetochore-microtubule attachments in cells and the successful distribution of the genetic material during cell division. The reversible inhibition of microtubule binding is an important aspect of the physiological error correction process. Small molecule inhibitors of protein-protein interactions involving Ndc80 are therefore highly desirable, both for mechanistic studies of chromosome segregation and also for their potential therapeutic value. Here, we report on a novel strategy to develop rationally designed inhibitors of the Ndc80 Calponin-homology domain using Supramolecular Chemistry. With a multiple-click approach, lysine-specific molecular tweezers were assembled to form covalently fused dimers to pentamers with a different overall size and preorganization/stiffness. We identified two dimers and a trimer as efficient Ndc80 CH-domain binders and have shown that they disrupt the interaction between Ndc80 and microtubules at low micromolar concentrations without affecting microtubule dynamics. NMR spectroscopy allowed us to identify the biologically important lysine residues 160 and 204 as preferred tweezer interaction sites. Enhanced sampling molecular dynamics simulations provided a rationale for the binding mode of multivalent tweezers and the role of pre-organization and secondary interactions in targeting multiple lysine residues across a protein surface.

Insights into substrate transport and water permeation in the mycobacterial transporter MmpL3

Mycobacteria, such as Mycobacterium tuberculosis, are characterized by a uniquely thick and waxy cell envelope that consists of two membranes, with a variety of mycolates comprising their outer membrane (OM). The protein Mycobacterial membrane protein Large 3 (MmpL3) is responsible for the transport of a primary OM component, trehalose monomycolate (TMM), from the inner (cytoplasmic) membrane (IM) to the periplasmic space, a process driven by the proton gradient. Although multiple structures of MmpL3 with bound substrates have been solved, the exact pathway(s) for TMM or proton transport remains elusive. Here, employing molecular dynamics simulations we investigate putative pathways for either transport species. We hypothesized that MmpL3 will cycle through similar conformational states as the related transporter AcrB, which we used as targets for modeling the conformation of MmpL3. A continuous water pathway through the transmembrane region was found in one of these states, illustrating a putative pathway for protons. Additional equilibrium simulations revealed that TMM can diffuse from the membrane into a binding pocket in MmpL3 spontaneously. We also found that acetylation of TMM, which is required for transport, makes it more stable within MmpL3's periplasmic cavity compared with the unacetylated form.

Unraveling the catalytic mechanism of SARS-CoV-2 papain-like protease with allosteric modulation of C270 mutation using multiscale computational approaches

Papain-like protease (PL^pro) is a promising therapeutic target against SARS-CoV-2, but its restricted S1/S2 subsites pose an obstacle in developing active site-directed inhibitors. We have recently identified C270 as a novel covalent allosteric site for SARS-CoV-2 PL^pro inhibitors. Here we present a theoretical investigation of the proteolysis reaction catalyzed by the wild-type SARS-CoV-2 PL^pro as well as the C270R mutant. Enhanced sampling MD simulations were first performed to explore the influence of C270R mutation on the protease dynamics, and sampled thermodynamically favorable conformations were then submitted to MM/PBSA and QM/MM MD simulations for thorough characterization of the protease-substrate binding and covalent reactions. The disclosed proteolysis mechanism of PL^pro, as characterized by the occurrence of proton transfer from the catalytic C111 to H272 prior to the substrate binding and with deacylation being the rate-determining step of the whole proteolysis process, is not completely identical to that of the 3C-like protease, another key cysteine protease of coronaviruses. The C270R mutation alters the structural dynamics of the BL2 loop that indirectly impairs the catalytic function of H272 and reduces the binding of the substrate with the protease, ultimately showing an inhibitory effect on PL^pro. Together, these results provide a comprehensive understanding at the atomic level of the key aspects of SARS-CoV-2 PL^pro proteolysis, including the catalytic activity allosterically regulated by C270 modification, which is crucial to the follow-up inhibitor design and development.

Deep Boosted Molecular Dynamics (DBMD): Accelerating molecular simulations with Gaussian boost potentials generated using probabilistic Bayesian deep neural network

We have developed a new Deep Boosted Molecular Dynamics (DBMD) method. Probabilistic Bayesian neural network models were implemented to construct boost potentials that exhibit Gaussian distribution with minimized anharmonicity, thereby allowing for accurate energetic reweighting and enhanced sampling of molecular simulations. DBMD was demonstrated on model systems of alanine dipeptide and the fast-folding protein and RNA structures. For alanine dipeptide, 30ns DBMD simulations captured up to 83-125 times more backbone dihedral transitions than 1μs conventional molecular dynamics (cMD) simulations and were able to accurately reproduce the original free energy profiles. Moreover, DBMD sampled multiple folding and unfolding events within 300ns simulations of the chignolin model protein and identified low-energy conformational states comparable to previous simulation findings. Finally, DBMD captured a general folding pathway of three hairpin RNAs with the GCAA, GAAA, and UUCG tetraloops. Based on Deep Learning neural network, DBMD provides a powerful and generally applicable approach to boosting biomolecular simulations. DBMD is available with open source in OpenMM at https://github.com/MiaoLab20/DBMD/.

Nanodisc-embedded cytochrome P450 P3A4 binds diverse ligands by distributing conformational dynamics to its flexible elements

Cytochrome P450 3A4 (CYP3A4) metabolizes a wide range of drugs and toxins. Interactions of CYP3A4 with ligands are difficult to predict due to promiscuity and conformational flexibility. To better understand CYP3A4 conformational responses to ligands we use hydrogen deuterium exchange mass spectrometry (HDX-MS) to investigate the effect of ligands on nanodisc-embedded CYP3A4. For a subset of CYP3A4-ligand complexes, differences in the low-frequency modes derived by principal component analyses of molecular dynamics trajectories mirrored the HDX-MS results. The effects of ligands are distributed to flexible elements of CYP3A4 between stretches of secondary structure. The largest effects occur in the F- and G-helices, where most ligands increase the flexibility of the F-helix and connecting loops and decrease the flexibility of the C-term of the G-helix. Most ligands affect the E-F-G, CD and HI regions of the protein. Ligand-dependent differences are observed in the A"-A' loop, BC region, E-helix, K-β1 region, proximal loop, and C-term loop. Correlated HDX responses were observed in the CD region and the C-term of the G-helix that were most pronounced for Type II ligands. Collectively, the HDX and molecular dynamics results suggest that CYP3A4 accommodates diverse binding partners by propagating local backbone fluctuations from the binding site onto the flexible regions of the enzyme via long-range interactions that are differentially modulated by ligands. In contrast to the paradigm wherein ligands decrease protein dynamics at their binding site, a wide range of ligands modestly increase CYP3A4 dynamics throughout the protein including effects remote from the active site.

Molecular dynamics study reveals key disruptors of MEIG1-PACRG interaction

Interactions between the meiosis-expressed gene 1 (MEIG1) and Parkin co-regulated gene (PACRG) protein are critical in the formation of mature sperm cells. Targeting either MEIG1 or PACRG protein could be a contraceptive strategy. The W50A and Y68A mutations on MEIG1 are known to interrupt the MEIG1-PACRG interactions resulting in defective sperm cells. However, the details about how the mutants disrupt the protein-protein binding are not clear. In this study, we reveal insights on MEIG1 and PACRG protein dynamics by applying Gaussian-accelerated molecular dynamics (GaMD) simulations and post-GaMD analysis. Our results show that the mutations destabilize the protein-protein interfacial interaction. The effect of the Y68A mutation is more significant than W50A as Y68 forms stronger polar interactions with PACRG. Because both human and mouse models demonstrate similar dynamic properties, the findings from mouse proteins can be applied to the human system. Moreover, we report a potential ligand binding pocket on the MEIG1 and PACRG interaction surface that could be a target for future drug design to inhibit the MEIG1-PACRG interaction. PACRG shows more qualified pockets along the protein-protein interface, implying that it is a better target than MEIG1. Our work provides a fundamental understanding of MEIG1 and PACRG protein dynamics, paving the way for drug discovery in male-based contraception.

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.

Elucidation of a dynamic interplay between a beta-2 adrenergic receptor, its agonist, and stimulatory G protein

G protein-coupled receptors (GPCRs) represent the largest group of membrane receptors for transmembrane signal transduction. Ligand-induced activation of GPCRs triggers G protein activation followed by various signaling cascades. Understanding the structural and energetic determinants of ligand binding to GPCRs and GPCRs to G proteins is crucial to the design of pharmacological treatments targeting specific conformations of these proteins to precisely control their signaling properties. In this study, we focused on interactions of a prototypical GPCR, beta-2 adrenergic receptor (β2AR), with its endogenous agonist, norepinephrine (NE), and the stimulatory G protein (Gs). Using molecular dynamics (MD) simulations, we demonstrated the stabilization of cationic NE, NE(+), binding to β2AR by Gs protein recruitment, in line with experimental observations. We also captured the partial dissociation of the ligand from β2AR and the conformational interconversions of Gs between closed and open conformations in the NE(+)-β2AR-Gs ternary complex while it is still bound to the receptor. The variation of NE(+) binding poses was found to alter Gs α subunit (Gsα) conformational transitions. Our simulations showed that the interdomain movement and the stacking of Gsα α1 and α5 helices are significant for increasing the distance between the Gsα and β2AR, which may indicate a partial dissociation of Gsα The distance increase commences when Gsα is predominantly in an open state and can be triggered by the intracellular loop 3 (ICL3) of β2AR interacting with Gsα, causing conformational changes of the α5 helix. Our results help explain molecular mechanisms of ligand and GPCR-mediated modulation of G protein activation.

Elucidating the molecular determinants in the process of gastrin C-terminal pentapeptide amide end activating cholecystokinin 2 receptor by Gaussian accelerated molecular dynamics simulations

Gastrin plays important role in stimulating the initiation and development of many gastrointestinal diseases through interacting with the cholecystokinin 2 receptor (CCK2R). The smallest bioactive unit of gastrin activating CCK2R is the C-terminal tetrapeptide capped with an indispensable amide end. Understanding the mechanism of this smallest bioactive unit interacting with CCK2R on a molecular basis could provide significant insights for designing CCK2R antagonists, which can be used to treat gastrin-related diseases. To this end, we performed extensive Gaussian accelerated molecular dynamics simulations to investigate the interaction between gastrin C-terminal pentapeptide capped with/without amide end and CCK2R. The amide cap influences the binding modes of the pentapeptide with CCK2R by weakening the electrostatic attractions between the C-terminus of the pentapeptide and basic residues near the extracellular domain in CCK2R. The C-terminus with the amide cap penetrates into the transmembrane domain of CCK2R while floating at the extracellular domain without the amide cap. Different binding modes induced different conformational dynamics of CCK2R. Residue pairs in CCK2R had stronger correlated motions when binding with the amidated pentapeptide. Key residues and interactions important for CCK2R binding with the amidated pentagastrin were also identified. Our results provide molecular insights into the determinants of the bioactive unit of gastrin activating CCK2R, which would be of great help for the design of CCK2R antagonists.

Allosteric Boost by TAB1 on the TAK1 Kinase Favorably Sculpts the Thermodynamic Landscape of Activation

The intricate mechanisms of allosteric regulation in kinases are of general interest to the scientific community for potential therapeutic implications. However, the diversity among kinases and their regulatory routes requires a case-by-case study to widen the repertoire of known mechanisms. The present study achieves this by understanding TAK1 kinase activation by TAB1 as a model phenomenon for the first time. Despite the known capacity of TAK1 to switch between its inactive ("DFG-out") and active-like ("DFG-in") conformations, the questionable role of TAB1 in offering an energetic favor to this has been addressed here using sequential combination of enhanced sampling methods like targeted molecular dynamics (TMD) and Gaussian accelerated molecular dynamics (GaMD). It reveals how a minimal domain of TAB1 sufficiently acts like a "catalytic gear" by favorably sculpting TAK1's thermodynamic landscape (potential of mean force in 2D) that accelerates "in"-"out" conformational switching of the conserved DFG motif. Standard molecular dynamics simulations (∼5 μs) reveal that TAB1 fascinatingly exploits the "lever-like" αF helix of TAK1 kinase domain to remotely propel the DFG motif via subtle helical "unfolding-folding" modifications within the kinase activation loop. The presence of two charged residues on terminal poles of αF helix imparts it, with this unique "lever-like" utility, and this turns out to be one important signature of co-evolution between TAK1 and TAB1. The entire mechanism of TAB1's impact transduction, which is found to be analogous to the moves in the popular "Chinese checker" game, gives a clear proof of the "dynamics-driven allostery" concept in kinases. The findings further benchmark TAK1's known autophosphorylation capacity. A novel insight into kinase allostery is thus provided, which potentiates investigation of similar capacities in other kinases.

Lipid membrane mimetics and oligomerization tune functional properties of proteorhodopsin

The functional properties of proteorhodopsin (PR) have been found to be strongly modulated by oligomeric distributions and lipid membrane mimetics. This study aims to distinguish and explain their effects by investigating how oligomer formation impacts PR's function of proton transport in lipid-based membrane mimetic environments. We find that PR forms stable hexamers and pentamers in both E. coli membranes and synthetic liposomes. Compared with the monomers, the photocycle kinetics of PR oligomers is ∼2 and ∼4.5 times slower for transitions between the K and M and the M and N photointermediates, respectively, indicating that oligomerization significantly slows PR's rate of proton transport in liposomes. In contrast, the apparent pKa of the key proton acceptor residue D97 (pKaD97) of liposome-embedded PR persists at 6.2-6.6, regardless of cross-protomer modulation of D97, suggesting that the liposome environment helps maintain PR's functional activity at neutral pH. By comparison, when extracted directly from E. coli membranes into styrene-maleic acid lipid particles, the pKaD97 of monomer-enriched E50Q PR drastically increases to 8.9, implying that there is a very low active PR population at neutral pH to engage in PR's photocycle. These findings demonstrate that oligomerization impacts PR's photocycle kinetics, while lipid-based membrane mimetics strongly affect PR's active population via different mechanisms.

CHARMM-GUI Enhanced Sampler for various collective variables and enhanced sampling methods

Enhanced sampling methodologies modifying underlying Hamiltonians can be used for the systems with a rugged potential energy surface that makes it hard to observe convergence using conventional unbiased molecular dynamics (MD) simulations. We present CHARMM-GUI Enhanced Sampler, a web-based tool to prepare various enhanced sampling simulations inputs with user-selected collective variables (CVs). Enhanced Sampler provides inputs for the following nine methods: accelerated MD, Gaussian accelerated MD, conformational flooding, metadynamics, adaptive biasing force, steered MD, temperature replica exchange MD, replica exchange solute tempering 2, and replica exchange umbrella sampling for the method-implemented MD packages including AMBER, CHARMM, GENESIS, GROMACS, NAMD, and OpenMM. Users only need to select a group of atoms via intuitive web-implementation in order to define commonly used nine CVs of interest: center of mass based distance, angle, dihedral, root-mean-square-distance, radius of gyration, distance projected on axis, two types of angles projected on axis, and coordination numbers. The enhanced sampling methods are tested with several biological systems to illustrate their efficiency over conventional MD. Enhanced Sampler with carefully optimized system-dependent parameters will help users to get meaningful results from their enhanced sampling simulations.

A genetically encoded BRET-based SARS-CoV-2 Mpro protease activity sensor

The main protease, Mpro, is critical for SARS-CoV-2 replication and an appealing target for designing anti-SARS-CoV-2 agents. Therefore, there is a demand for the development of improved sensors to monitor its activity. Here, we report a pair of genetically encoded, bioluminescence resonance energy transfer (BRET)-based sensors for detecting Mpro proteolytic activity in live cells as well as in vitro. The sensors were generated by sandwiching peptides containing the Mpro N-terminal autocleavage sites, either AVLQSGFR (short) or KTSAVLQSGFRKME (long), in between the mNeonGreen and NanoLuc proteins. Co-expression of the sensors with Mpro in live cells resulted in their cleavage while mutation of the critical C145 residue (C145A) in Mpro completely abrogated their cleavage. Additionally, the sensors recapitulated the inhibition of Mpro by the well-characterized pharmacological agent GC376. Further, in vitro assays with the BRET-based Mpro sensors revealed a molecular crowding-mediated increase in the rate of Mpro activity and a decrease in the inhibitory potential of GC376. The sensors developed here will find direct utility in studies related to drug discovery targeting the SARS-CoV-2 Mpro and functional genomics application to determine the effect of sequence variation in Mpro.

Advanced Molecular Tweezers with Lipid Anchors against SARS-CoV-2 and Other Respiratory Viruses

The COVID-19 pandemic caused by SARS-CoV-2 presents a global health emergency. Therapeutic options against SARS-CoV-2 are still very limited but urgently required. Molecular tweezers are supramolecular agents that destabilize the envelope of viruses resulting in a loss of viral infectivity. Here, we show that first-generation tweezers, CLR01 and CLR05, disrupt the SARS-CoV-2 envelope and abrogate viral infectivity. To increase the antiviral activity, a series of 34 advanced molecular tweezers were synthesized by insertion of aliphatic or aromatic ester groups on the phosphate moieties of the parent molecule CLR01. A structure-activity relationship study enabled the identification of tweezers with a markedly enhanced ability to destroy lipid bilayers and to suppress SARS-CoV-2 infection. Selected tweezer derivatives retain activity in airway mucus and inactivate the SARS-CoV-2 wildtype and variants of concern as well as respiratory syncytial, influenza, and measles viruses. Moreover, inhibitory activity of advanced tweezers against respiratory syncytial virus and SARS-CoV-2 was confirmed in mice. Thus, potentiated tweezers are broad-spectrum antiviral agents with great prospects for clinical development to combat highly pathogenic viruses.

Gaussian Accelerated Molecular Dynamics in OpenMM

Gaussian accelerated molecular dynamics (GaMD) is a computational technique that provides both unconstrained enhanced sampling and free energy calculations of biomolecules. Here, we present the implementation of GaMD in the OpenMM simulation package and validate it on model systems of alanine dipeptide and RNA folding. For alanine dipeptide, 30 ns GaMD production simulations reproduced free energy profiles of 1000 ns conventional molecular dynamics (cMD) simulations. In addition, GaMD simulations captured the folding pathways of three hyperstable RNA tetraloops (UUCG, GCAA, and CUUG) and binding of the rbt203 ligand to the HIV-1 Tar RNA, both of which involved critical electrostatic interactions such as hydrogen bonding and base stacking. Together with previous implementations, GaMD in OpenMM will allow for wider applications in simulations of proteins, RNA, and other biomolecules.

Potential Use of Tea Tree Oil as a Disinfectant Agent against Coronaviruses: A Combined Experimental and Simulation Study

The COVID-19 pandemic has highlighted the relevance of proper disinfection procedures and renewed interest in developing novel disinfectant materials as a preventive strategy to limit SARS-CoV-2 contamination. Given its widely known antibacterial, antifungal, and antiviral properties, Melaleuca alternifolia essential oil, also named Tea tree oil (TTO), is recognized as a potential effective and safe natural disinfectant agent. In particular, the proposed antiviral activity of TTO involves the inhibition of viral entry and fusion, interfering with the structural dynamics of the membrane and with the protein envelope components. In this study, for the first time, we demonstrated the virucidal effects of TTO against the feline coronavirus (FCoVII) and the human coronavirus OC43 (HCoV-OC43), both used as surrogate models for SARS-CoV-2. Then, to atomistically uncover the possible effects exerted by TTO compounds on the outer surface of the SARS-CoV-2 virion, we performed Gaussian accelerated Molecular Dynamics simulations of a SARS-CoV-2 envelope portion, including a complete model of the Spike glycoprotein in the absence or presence of the three main TTO compounds (terpinen-4-ol, γ-terpinene, and 1,8-cineole). The obtained results allowed us to hypothesize the mechanism of action of TTO and its possible use as an anti-coronavirus disinfectant agent.

Mechanism of Peptide Agonist Binding in CXCR4 Chemokine Receptor

Chemokine receptors are key G-protein-coupled receptors (GPCRs) that control cell migration in immune system responses, development of cardiovascular and central nervous systems, and numerous diseases. In particular, the CXCR4 chemokine receptor promotes metastasis, tumor growth and angiogenesis in cancers. CXCR4 is also used as one of the two co-receptors for T-tropic HIV-1 entry into host cells. Therefore, CXCR4 serves as an important therapeutic target for treating cancers and HIV infection. Apart from the CXCL12 endogenous peptide agonist, previous studies suggested that the first 17 amino acids of CXCL12 are sufficient to activate CXCR4. Two 17-residue peptides with positions 1-4 mutated to RSVM and ASLW functioned as super and partial agonists of CXCR4, respectively. However, the mechanism of peptide agonist binding in CXCR4 remains unclear. Here, we have investigated this mechanism through all-atom simulations using a novel Peptide Gaussian accelerated molecular dynamics (Pep-GaMD) method. The Pep-GaMD simulations have allowed us to explore representative binding conformations of each peptide and identify critical low-energy states of CXCR4 activated by the super versus partial peptide agonists. Our simulations have provided important mechanistic insights into peptide agonist binding in CXCR4, which are expected to facilitate rational design of new peptide modulators of CXCR4 and other chemokine receptors.

Protein-Protein Interaction-Gaussian Accelerated Molecular Dynamics (PPI-GaMD): Characterization of Protein Binding Thermodynamics and Kinetics

Protein-protein interactions (PPIs) play key roles in many fundamental biological processes such as cellular signaling and immune responses. However, it has proven challenging to simulate repetitive protein association and dissociation in order to calculate binding free energies and kinetics of PPIs due to long biological timescales and complex protein dynamics. To address this challenge, we have developed a new computational approach to all-atom simulations of PPIs based on a robust Gaussian accelerated molecular dynamics (GaMD) technique. The method, termed "PPI-GaMD", selectively boosts interaction potential energy between protein partners to facilitate their slow dissociation. Meanwhile, another boost potential is applied to the remaining potential energy of the entire system to effectively model the protein's flexibility and rebinding. PPI-GaMD has been demonstrated on a model system of the ribonuclease barnase interactions with its inhibitor barstar. Six independent 2 μs PPI-GaMD simulations have captured repetitive barstar dissociation and rebinding events, which enable calculations of the protein binding thermodynamics and kinetics simultaneously. The calculated binding free energies and kinetic rate constants agree well with the experimental data. Furthermore, PPI-GaMD simulations have provided mechanistic insights into barstar binding to barnase, which involves long-range electrostatic interactions and multiple binding pathways, being consistent with previous experimental and computational findings of this model system. In summary, PPI-GaMD provides a highly efficient and easy-to-use approach for binding free energy and kinetics calculations of PPIs.

GLOW: A Workflow Integrating Gaussian-Accelerated Molecular Dynamics and Deep Learning for Free Energy Profiling

We introduce a Gaussian-accelerated molecular dynamics (GaMD), deep learning (DL), and free energy profiling workflow (GLOW) to predict molecular determinants and map free energy landscapes of biomolecules. All-atom GaMD-enhanced sampling simulations are first performed on biomolecules of interest. Structural contact maps are then calculated from GaMD simulation frames and transformed into images for building DL models using a convolutional neural network. Important structural contacts are further determined from DL models of attention maps of the structural contact gradients, which allow us to identify the system reaction coordinates. Finally, free energy profiles are calculated for the selected reaction coordinates through energetic reweighting of the GaMD simulations. We have also successfully demonstrated GLOW for the characterization of activation and allosteric modulation of a G protein-coupled receptor, using the adenosine A₁ receptor (A₁AR) as a model system. GLOW findings are highly consistent with previous experimental and computational studies of the A₁AR, while also providing further mechanistic insights into the receptor function. In summary, GLOW provides a systematic approach to mapping free energy landscapes of biomolecules. The GLOW workflow and its user manual can be downloaded at http://miaolab.org/GLOW.

Structure of the 5' untranslated region in SARS-CoV-2 genome and its specific recognition by innate immune system via the human oligoadenylate synthase 1

2'-5'-Oligoadenylate synthetase 1 (OAS1) is one of the key enzymes driving the innate immune system response to SARS-CoV-2 infection whose activity has been related to COVID-19 severity. OAS1 is a sensor of endogenous RNA that triggers the 2'-5'-oligoadenylate/RNase L pathway. Upon SARS-CoV-2 infection, OAS1 is responsible for the recognition of viral RNA and has been shown to possess a particularly high sensitivity for the 5'-untranslated (5'-UTR) RNA region, which is organized in a double-strand stem loop motif (SL1). Here we report the structure of the SL1/OAS1 complex also rationalizing the high affinity for OAS1.

The role of DNA nanostructures in the catalytic properties of an allosterically regulated protease

DNA-scaffolded enzymes typically show altered kinetic properties; however, the mechanism behind this phenomenon is still poorly understood. We address this question using thrombin, a model of allosterically regulated serine proteases, encaged into DNA origami cavities with distinct structural and electrostatic features. We compare the hydrolysis of substrates that differ only in their net charge due to a terminal residue far from the cleavage site and presumably involved in the allosteric activation of thrombin. Our data show that the reaction rate is affected by DNA/substrate electrostatic interactions, proportionally to the degree of DNA/enzyme tethering. For substrates of opposite net charge, this leads to an inversion of the catalytic response of the DNA-scaffolded thrombin when compared to its freely diffusing counterpart. Hence, by altering the electrostatic environment nearby the encaged enzyme, DNA nanostructures interfere with charge-dependent mechanisms of enzyme-substrate recognition and may offer an alternative tool to regulate allosteric processes through spatial confinement.

Binding Analysis Using Accelerated Molecular Dynamics Simulations and Future Perspectives

Biomolecular recognition such as binding of small molecules, nucleic acids, peptides and proteins to their target receptors plays key roles in cellular function and has been targeted for therapeutic drug design. Molecular dynamics (MD) is a computational approach to analyze these binding processes at an atomistic level, which provides valuable understandings of the mechanisms of biomolecular recognition. However, the rather slow biomolecular binding events often present challenges for conventional MD (cMD), due to limited simulation timescales (typically over hundreds of nanoseconds to tens of microseconds). In this regard, enhanced sampling methods, particularly accelerated MD (aMD), have proven useful to bridge the gap and enable all-atom simulations of biomolecular binding events. Here, we will review the recent method developments of Gaussian aMD (GaMD), ligand GaMD (LiGaMD) and peptide GaMD (Pep-GaMD), which have greatly expanded our capabilities to simulate biomolecular binding processes. Spontaneous binding of various biomolecules to their receptors has been successfully simulated by GaMD. Microsecond LiGaMD and Pep-GaMD simulations have captured repetitive binding and dissociation of small-molecule ligands and highly flexible peptides, and thus enabled ligand/peptide binding thermodynamics and kinetics calculations. We will also present relevant application studies in simulations of important drug targets and future perspectives for rational computer-aided drug design.

Challenges and frontiers of computational modelling of biomolecular recognition

Biomolecular recognition including binding of small molecules, peptides and proteins to their target receptors plays a key role in cellular function and has been targeted for therapeutic drug design. However, the high flexibility of biomolecules and slow binding and dissociation processes have presented challenges for computational modelling. Here, we review the challenges and computational approaches developed to characterise biomolecular binding, including molecular docking, molecular dynamics simulations (especially enhanced sampling) and machine learning. Further improvements are still needed in order to accurately and efficiently characterise binding structures, mechanisms, thermodynamics and kinetics of biomolecules in the future.

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

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

Mechanism of RNA recognition by a Musashi RNA-binding protein

The Musashi RNA-binding proteins (RBPs) regulate translation of target mRNAs and maintenance of cell stemness and tumorigenesis. Musashi-1 (MSI1), long considered as an intestinal and neural stem cell marker, has been more recently found to be over expressed in many cancers. It has served as an important drug target for treating acute myeloid leukemia and solid tumors such as ovarian, colorectal and bladder cancer. One of the reported binding targets of MSI1 is Numb, a negative regulator of the Notch signaling. However, the dynamic mechanism of Numb RNA binding to MSI1 remains unknown, largely hindering effective drug design targeting this critical interaction. Here, we have performed extensive all-atom microsecond-timescale simulations using a robust Gaussian accelerated molecular dynamics (GaMD) method, which successfully captured multiple times of spontaneous and highly accurate binding of the Numb RNA from bulk solvent to the MSI1 protein target site. GaMD simulations revealed that Numb RNA binding to MSI1 involved largely induced fit in both the RNA and protein. The simulations also identified important low-energy intermediate conformational states during RNA binding, in which Numb interacted mainly with the β2-β3 loop and C terminus of MSI1. The mechanistic understanding of RNA binding obtained from our GaMD simulations is expected to facilitate rational structure-based drug design targeting MSI1 and other RBPs.

Molecular Insights into Site-Specific Interferon-α2a Bioconjugates Originated from PEG, LPG, and PEtOx

Conjugation of biologics with polymers modulates their pharmacokinetics, with polyethylene glycol (PEG) as the gold standard. We compared alternative polymers and two types of cyclooctyne linkers (BCN/DBCO) for bioconjugation of interferon-α2a (IFN-α2a) using 10 kDa polymers including linear mPEG, poly(2-ethyl-2-oxazoline) (PEtOx), and linear polyglycerol (LPG). IFN-α2a was azide functionalized via amber codon expansion and bioorthogonally conjugated to all cyclooctyne linked polymers. Polymer conjugation did not impact IFN-α2a's secondary structure and only marginally reduced IFN-α2a's bioactivity. In comparison to PEtOx, the LPG polymer attached via the less rigid cyclooctyne linker BCN was found to stabilize IFN-α2a against thermal stress. These findings were further detailed by molecular modeling studies which showed a modulation of protein flexibility upon PEtOx conjugation and a reduced amount of protein native contacts as compared to PEG and LPG originated bioconjugates. Polymer interactions with IFN-α2a were further assessed via a limited proteolysis (LIP) assay, which resulted in comparable proteolytic cleavage patterns suggesting weak interactions with the protein's surface. In conclusion, both PEtOx and LPG bioconjugates resulted in a similar biological outcome and may become promising PEG alternatives for bioconjugation.

Molecular dynamics study of enhanced autophosphorylation by S904F mutation of the RET kinase domain

The aberrant kinase activity of RET (rearranged during transfection), a transmembrane tyrosine kinase, is associated with human cancer. A point mutation caused by the replacement of solvent-front hydrophilic S904, located on the activation loop (A-loop), with a bulky hydrophobic phenylalanine residue can induce resistance to the type I kinase inhibitor vandetanib. A possible mechanism of this drug resistance is the release of a cis-autoinhibited conformation of RET for autophosphorylation, which activates RET kinase. Because the association between S904F mutation and enhanced autophosphorylation is unclear, we conducted molecular modeling analysis to compare unphosphorylated apo wild-type and S904F mutant structures. The structural compactness of the A-loop promoted ATP binding. When the A-loop is extended, the αC helix moves toward the glycine-rich loop, resulting in the protrusion of F735. The extruded F735 connects with E734 and R912 and constrains the ATP pocket entrance. Contrarily, a contracted A-loop pulls the αC helix away from the glycine-rich loop, burying F734 and making the ATP pocket accessible. The mutated F904 stabilizes the contracted A-loop and releases the autoinhibited conformation of RET, thereby facilitating autophosphorylation. We also simulated two ATP-bound systems. The binding free energies of ATP, estimated through the molecular mechanics with a generalized Born and surface area solvation approach, revealed that the S904F mutant was bound more tightly than was the wild type with the ATP. The findings support the premise of autophosphorylation promotion in the S904F mutant.

Multiscale computational study of ligand binding pathways: Case of p38 MAP kinase and its inhibitors

Protein kinases are one of the most important drug targets in the past 10 years. Understanding the inhibitor association processes will profoundly impact new binder designs with preferred binding kinetics. However, after more than a decade of effort, a complete atomistic-level study of kinase inhibitor binding pathways is still lacking. As all kinases share a similar scaffold, we used p38 kinase as a model system to investigate the conformational dynamics and free energy transition of inhibitor binding toward kinases. Two major kinase conformations, Asp-Phe-Gly (DFG)-in and DFG-out, and three types of inhibitors, type I, II, and III, were thoroughly investigated in this work. We performed Brownian dynamics simulations and up to 340 μs Gaussian-accelerated molecular dynamics simulations to capture the inhibitor binding paths and a series of conformational transitions of the p38 kinase from its apo to inhibitor-bound form. Eighteen successful binding trajectories, including all types of inhibitors, are reported herein. Our simulations suggest a mechanism of inhibitor recruitment, a faster ligand association step to a pre-existing DFG-in/DFG-out p38 protein, followed by a slower molecular rearrangement step to adjust the protein-ligand conformation followed by a shift in the energy landscape to reach the final bound state. The ligand association processes also reflect the energetic favor of type I and type II/III inhibitor binding through ATP and allosteric channels, respectively. These different binding routes are directly responsible for the fast (type I binders) and slow (type II/III binders) kinetics of different types of p38 inhibitors. Our findings also echo the recent study of p38 inhibitor dissociation, implying that ligand unbinding could undergo a reverse path of binding, and both processes share similar metastates. This study deepens the understanding of molecular and energetic features of kinase inhibitor-binding processes and will inspire future drug development from a kinetic point of view.

Specific Engineered G Protein Coupling to Histamine Receptors Revealed from Cellular Assay Experiments and Accelerated Molecular Dynamics Simulations

G protein-coupled receptors (GPCRs) are targets of extracellular stimuli and hence occupy a key position in drug discovery. By specific and not yet fully elucidated coupling profiles with α subunits of distinct G protein families, they regulate cellular responses. The histamine H2 and H4 receptors (H2R and H4R) are prominent members of Gs- and Gi-coupled GPCRs. Nevertheless, promiscuous G protein and selective Gi signaling have been reported for the H2R and H4R, respectively, the molecular mechanism of which remained unclear. Using a combination of cellular experimental assays and Gaussian accelerated molecular dynamics (GaMD) simulations, we investigated the coupling profiles of the H2R and H4R to engineered mini-G proteins (mG). We obtained coupling profiles of the mGs, mGsi, or mGsq proteins to the H2R and H4R from the mini-G protein recruitment assays using HEK293T cells. Compared to H2R-mGs expressing cells, histamine responses were weaker (pEC50, Emax) for H2R-mGsi and -mGsq. By contrast, the H4R selectively bound to mGsi. Similarly, in all-atom GaMD simulations, we observed a preferential binding of H2R to mGs and H4R to mGsi revealed by the structural flexibility and free energy landscapes of the complexes. Although the mG α5 helices were consistently located within the HR binding cavity, alternative binding orientations were detected in the complexes. Due to the specific residue interactions, all mG α5 helices of the H2R complexes adopted the Gs-like orientation toward the receptor transmembrane (TM) 6 domain, whereas in H4R complexes, only mGsi was in the Gi-like orientation toward TM2, which was in agreement with Gs- and Gi-coupled GPCRs structures resolved by X-ray/cryo-EM. These cellular and molecular insights support (patho)physiological profiles of the histamine receptors, especially the hitherto little studied H2R function in the brain, as well as of the pharmacological potential of H4R selective drugs.

Gaussian accelerated molecular dynamics (GaMD): principles and applications

Gaussian accelerated molecular dynamics (GaMD) is a robust computational method for simultaneous unconstrained enhanced sampling and free energy calculations of biomolecules. It works by adding a harmonic boost potential to smooth biomolecular potential energy surface and reduce energy barriers. GaMD greatly accelerates biomolecular simulations by orders of magnitude. Without the need to set predefined reaction coordinates or collective variables, GaMD provides unconstrained enhanced sampling and is advantageous for simulating complex biological processes. The GaMD boost potential exhibits a Gaussian distribution, thereby allowing for energetic reweighting via cumulant expansion to the second order (i.e., "Gaussian approximation"). This leads to accurate reconstruction of free energy landscapes of biomolecules. Hybrid schemes with other enhanced sampling methods, such as the replica exchange GaMD (rex-GaMD) and replica exchange umbrella sampling GaMD (GaREUS), have also been introduced, further improving sampling and free energy calculations. Recently, new "selective GaMD" algorithms including the ligand GaMD (LiGaMD) and peptide GaMD (Pep-GaMD) enabled microsecond simulations to capture repetitive dissociation and binding of small-molecule ligands and highly flexible peptides. The simulations then allowed highly efficient quantitative characterization of the ligand/peptide binding thermodynamics and kinetics. Taken together, GaMD and its innovative variants are applicable to simulate a wide variety of biomolecular dynamics, including protein folding, conformational changes and allostery, ligand binding, peptide binding, protein-protein/nucleic acid/carbohydrate interactions, and carbohydrate/nucleic acid interactions. In this review, we present principles of the GaMD algorithms and recent applications in biomolecular simulations and drug design.

Anisotropy in mechanical unfolding of protein upon partner-assisted pulling and handle-assisted pulling

Proteins as force-sensors respond to mechanical cues and regulate signaling in physiology. Proteins commonly connect the source and response points of mechanical cues in two conformations, independent proteins in end-to-end geometry and protein complexes in handshake geometry. The force-responsive property of independent proteins in end-to-end geometry is studied extensively using single-molecule force spectroscopy (SMFS). The physiological significance of the complex conformations in force-sensing is often disregarded as mere surge protectors. However, with the potential of force-steering, protein complexes possess a distinct mechano-responsive property over individual force-sensors. To decipher, we choose a force-sensing protein, cadherin-23, from tip-link complex and perform SMFS using end-to-end geometry and handshake complex geometry. We measure higher force-resilience of cadherin-23 with preferential shorter extensions in handshake mode of pulling over the direct mode. The handshake geometry drives the force-response of cadherin-23 through different potential-energy landscapes than direct pulling. Analysis of the dynamic network structure of cadherin-23 under tension indicates narrow force-distributions among residues in cadherin-23 in direct pulling, resulting in low force-dissipation paths and low resilience to force. Overall, the distinct and superior mechanical responses of cadherin-23 in handshake geometry than single protein geometry highlight a probable evolutionary drive of protein-protein complexes as force-conveyors over independent ones.

Elucidation of Cryptic and Allosteric Pockets within the SARS-CoV-2 Main Protease

The SARS-CoV-2 pandemic has rapidly spread across the globe, posing an urgent health concern. Many quests to computationally identify treatments against the virus rely on in silico small molecule docking to experimentally determined structures of viral proteins. One limit to these approaches is that protein dynamics are often unaccounted for, leading to overlooking transient, druggable conformational states. Using Gaussian accelerated molecular dynamics to enhance sampling of conformational space, we identified cryptic pockets within the SARS-CoV-2 main protease, including some within regions far from the active site. These simulations sampled comparable dynamics and pocket volumes to conventional brute force simulations carried out on two orders of magnitude greater timescales.