Improved Reweighting of Accelerated Molecular Dynamics Simulations for Free Energy Calculation

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.

Combining Transition Path Sampling with Data-Driven Collective Variables through a Reactivity-Biased Shooting Algorithm

Rare event sampling is a central problem in modern computational chemistry research. Among the existing methods, transition path sampling (TPS) can generate unbiased representations of reaction processes. However, its efficiency depends on the ability to generate reactive trial paths, which in turn depends on the quality of the shooting algorithm used. We propose a new algorithm based on the shooting success rate, i.e., reactivity, measured as a function of a reduced set of collective variables (CVs). These variables are extracted with a machine learning approach directly from TPS simulations, using a multitask objective function. Iteratively, this workflow significantly improves the shooting efficiency without any prior knowledge of the process. In addition, the optimized CVs can be used with biased enhanced sampling methodologies to accurately reconstruct the free energy profiles. We tested the method on three different systems: a two-dimensional toy model, conformational transitions of alanine dipeptide, and hydrolysis of acetyl chloride in bulk water. In the latter, we integrated our workflow with an active learning scheme to learn a reactive machine learning-based potential, which allowed us to study the mechanism and free energy profile with an ab initio-like accuracy.

Probability Density Reweighting of High-Temperature Molecular Dynamics

Molecular dynamics (MD) simulation is a popular method for elucidating the structures and functions of biomolecules. However, exploring the conformational space, especially for large systems with slow transitions, often requires enhanced sampling methods. Although conducting MD at high temperatures provides a straightforward approach, resulting conformational ensembles diverge significantly from those at low temperatures. To address this discrepancy, we propose a novel probability density-based reweighting (PDR) method. PDR exhibits robust performance across four distinct systems, including a miniprotein, a cyclic peptide, a protein loop, and a protein-peptide complex. It accurately restores the conformational distributions at high temperatures to those at low temperatures. Additionally, we apply PDR to reweight previously studied high-T MD simulations of 12 protein-peptide complexes, enabling a comprehensive investigation of the conformational space of protein-peptide complexes.

Molecular Mechanism of Phosphorylation-Mediated Impacts on the Conformation Dynamics of GTP-Bound KRAS Probed by GaMD Trajectory-Based Deep Learning

The phosphorylation of different sites produces a significant effect on the conformational dynamics of KRAS. Gaussian accelerated molecular dynamics (GaMD) simulations were combined with deep learning (DL) to explore the molecular mechanism of the phosphorylation-mediated effect on conformational dynamics of the GTP-bound KRAS. The DL finds that the switch domains are involved in obvious differences in conformation contacts and suggests that the switch domains play a key role in the function of KRAS. The analyses of free energy landscapes (FELs) reveal that the phosphorylation of pY32, pY64, and pY137 leads to more disordered states of the switch domains than the wild-type (WT) KRAS and induces conformational transformations between the closed and open states. The results from principal component analysis (PCA) indicate that principal motions PC1 and PC2 are responsible for the closed and open states of the phosphorylated KRAS. Interaction networks were analyzed and the results verify that the phosphorylation alters interactions of GTP and magnesium ion Mg2+ with the switch domains. It is concluded that the phosphorylation pY32, pY64, and pY137 tune the activity of KRAS through changing conformational dynamics and interactions of the switch domains. We anticipated that this work could provide theoretical aids for deeply understanding the function of KRAS.

Ligand Gaussian accelerated Molecular Dynamics 3 (LiGaMD3): Improved Calculations of Binding Thermodynamics and Kinetics of Both Small Molecules and Flexible Peptides

Binding thermodynamics and kinetics play critical roles in drug design. However, it has proven challenging to efficiently predict ligand binding thermodynamics and kinetics of small molecules and flexible peptides using conventional Molecular Dynamics (cMD), due to limited simulation timescales. Based on our previously developed Ligand Gaussian accelerated Molecular Dynamics (LiGaMD) method, we present a new approach, termed "LiGaMD3", in which we introduce triple boosts into three individual energy terms that play important roles in small-molecule/peptide dissociation, rebinding and system conformational changes to improve the sampling efficiency of small-molecule/peptide interactions with target proteins. To validate the performance of LiGaMD3, MDM2 bound by a small molecule (Nutlin 3) and two highly flexible peptides (PMI and P53) were chosen as model systems. LiGaMD3 could efficiently capture repetitive small-molecule/peptide dissociation and binding events within 2 microsecond simulations. The predicted binding kinetic constant rates and free energies from LiGaMD3 agreed with available experimental values and previous simulation results. Therefore, LiGaMD3 provides a more general and efficient approach to capture dissociation and binding of both small-molecule ligand and flexible peptides, allowing for accurate prediction of their binding thermodynamics and kinetics.

The rotamer of the second-sphere histidine in AA9 lytic polysaccharide monooxygenase is pH dependent

Lytic polysaccharide monooxygenases (LPMOs) catalyze a reaction that is crucial for the biological decomposition of various biopolymers and for the industrial conversion of plant biomass. Despite the importance of LPMOs, the exact molecular-level nature of the reaction mechanism is still debated today. Here, we investigated the pH-dependent conformation of a second-sphere histidine (His) that we call the stacking histidine, which is conserved in fungal AA9 LPMOs and is speculated to assist catalysis in several of the LPMO reaction pathways. Using constant-pH and accelerated molecular dynamics simulations, we monitored the dynamics of the stacking His in different protonation states for both the resting Cu(II) and active Cu(I) forms of two fungal LPMOs. Consistent with experimental crystallographic and neutron diffraction data, our calculations suggest that the side chain of the protonated and positively charged form is rotated out of the active site toward the solvent. Importantly, only one of the possible neutral states of histidine (HIE state) is observed in the stacking orientation at neutral pH or when bound to cellulose. Our data predict that, in solution, the stacking His may act as a stabilizer (via hydrogen bonding) of the Cu(II)-superoxo complex after the LPMO-Cu(I) has reacted with O2 in solution, which, in fine, leads to H2O2 formation. Also, our data indicate that the HIE-stacking His is a poor acid/base catalyst when bound to the substrate and, in agreement with the literature, may play an important stabilizing role (via hydrogen bonding) during the peroxygenase catalysis. Our study reveals the pH titration midpoint values of the pH-dependent orientation of the stacking His should be considered when modeling and interpreting LPMO reactions, whether it be for classical LPMO kinetics or in industry-oriented enzymatic cocktails, and for understanding LPMO behavior in slightly acidic natural processes such as fungal wood decay.

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.

Revealing the key structural features promoting the helical conformation in algal polysaccharide carrageenan in solution

Carrageenans are industrially important polysaccharides with tunable viscoelastic and gelation properties. The function of polysaccharide depends on its conformation and chemical composition. However, the solution conformations of carrageenans are highly debated, and the structure-function relationship remains elusive. Here, we have studied the intrinsic conformational behavior of a series of carrageenan hexamers in solution, using extensive all-atom classical MD and enhanced sampling. Our findings comprehensively delineate that carrageenans containing the 3,6-anhydrous bridge (κ-C, ι-C, θ-C, and non-sulfated β-C) adopt compact helical structures as their predominant conformation in solution, whereas carrageenans without the bridge (μ-C, ν-C, and λ-C) remain as extended loosely packed helices; opposing the 'coil-to-helix' paradigm. Glycosidic linkages access a few allowed orientations. We hypothesize that the 3,6-anhydrous bridge, irrespective of carrageenan's sulfation pattern, is essential for stabilizing the helical conformation at the single-chain level. It provides necessary flexibility to the glycosidic linkage to sample conformations close to the experimentally derived helical structure and also prevents the sugar ring flipping. Sulfate groups mainly modify the chain stiffness due to steric and stereo-electronic effects and participate in hydrogen bonding. Such atomistic insights will be helpful for understanding the differential gelation mechanisms of carrageenans and fine-tuning polysaccharide backbone for various industrial applications.

Computational insights into the binding modes, keto-enol tautomerization and stereo-electronically controlled decarboxylation of oxaloacetate in the active site of macrophomate synthase

Oxaloacetic acid (OAA) is a β-ketocarboxylic acid, which plays an important role as an intermediate in some metabolic pathways, including the tricarboxylic acid cycle, gluconeogenesis and fatty acid biosynthesis. Animal studies have indicated that supplementing oxaloacetic acid shows an increase of lifespan and other substantial health benefits including mitochondrial DNA protection, and protection of retinal, neural and pancreatic tissues. Most of the chemical transformations of OAA in the metabolic pathways have been extensively studied; however, the understanding of decarboxylation of OAA at the atomic level is relatively lacking. Here, we carried out MD simulations and combined quantum mechanical/molecular mechanical (QM/MM) calculations as an example to systematically elucidate the binding modes, keto-enol tautomerization and decarboxylation of OAA in the active site of macrophomate synthase (MPS), which is a Mg(II)-dependent bifunctional enzyme that catalyzes both the decarboxylation of OAA and [4+2] cycloaddition of 2-pyrone with the decarboxylated intermediate of OAA (pyruvate enolate). On the basis of our calculations, it was found that the Mg2+-coordinated oxaloacetate may exist in enol forms and keto forms. The four keto forms can be transformed into each other by simply rotating the C2-C3 single bond, nevertheless, the keto-enol tautomerization strictly requires the assistance of pocket water molecules. In addition, the decarboxylation is stereo-electronically controlled, i.e., it is the relative orientation of the terminal carboxyl anion that determines the rate of decarboxylation. As such, the chemistry of oxaloacetate in the active site of MPS is complex. On one hand, the most stable binding mode (K-I) may undergo enol-keto tautomerization to isomerize to the enol form, which may further react with the second substrate; on the other hand, K-I may isomerize to another binding mode K-II to proceed decarboxylation to generate pyruvate enolate and CO2. Starting from K-I, the enol-keto tautomerization corresponds to a barrier of 16.2 kcal mol-1, whereas the decarboxylation is associated with an overall barrier of 19.7 kcal mol-1. These findings may provide useful information for understanding the chemistry of OAA and the catalysis of related enzymes, and they are basically in agreement with the available experimental kinetic data.

Studying the Collective Functional Response of a Receptor in Alchemical Ligand Binding Free Energy Simulations with Accelerated Solvation Layer Dynamics

Ligand binding free energy simulations (LB-FES) that involve sampling of protein functional conformations have been longstanding challenges in research on molecular recognition. Particularly, modeling of the conformational transition pathway and design of the heuristic biasing mechanism are severe bottlenecks for the existing enhanced configurational sampling (ECS) methods. Inspired by the key role of hydration in regulating conformational dynamics of macromolecules, this report proposes a novel ECS approach that facilitates binding-associated structural dynamics by accelerated hydration transitions in combination with the λ-exchange of free energy perturbation (FEP). Two challenging protein-ligand binding processes involving large configurational transitions of the receptor are studied, with hydration transitions at binding sites accelerated by Hamiltonian-simulated annealing of the hydration layer. Without the need for pathway analysis or ad hoc barrier flattening potential, LB-FES were performed with FEP/λ-exchange molecular dynamics simulation at a minor overhead for annealing of the hydration layer. The LB-FES studies showed that the accelerated rehydration significantly enhances the collective conformational transitions of the receptor, and convergence of binding affinity calculations is obtained at a sweet-spot simulation time scale. Alchemical LB-FES with the proposed ECS strategy is free from the effort of trial and error for the setup and realizes efficient on-the-fly sampling for the collective functional response of the receptor and bound water and therefore presents a practical approach to high-throughput screening in drug discovery.

Structure and function of an intermediate GPCR-Gαβγ complex

Unraveling the signaling roles of intermediate complexes is pivotal for G protein-coupled receptor (GPCR) drug development. Despite hundreds of GPCR-Gαβγ structures, these snapshots primarily capture the fully activated end-state complex. Consequently, a comprehensive understanding of the conformational transitions during GPCR activation and the roles of intermediate GPCR-G protein complexes in signaling remain elusive. Guided by a conformational landscape profiled by 19 F quantitative NMR ( 19 F-qNMR) and Molecular Dynamics (MD) simulations, we resolved the structure of an unliganded GPCR-G protein intermediate complex by blocking its transition to the fully activated end-state complex. More importantly, we presented direct evidence that the intermediate GPCR-Gαsβγ complex initiates a rate-limited nucleotide exchange without progressing to the fully activated end-state complex, thereby bridging a significant gap in our understanding the complexity of GPCR signaling. Understanding the roles of individual conformational states and their complexes in signaling efficacy and bias will help us to design drugs that discriminately target a disease-related conformation.

Mechanism of antagonist ligand binding to REV-ERBα

REV-ERBα, a therapeutically promising nuclear hormone receptor, plays a crucial role in regulating various physiological processes such as the circadian clock, inflammation, and metabolism. However, the availability of chemical probes to investigate the pharmacology of this receptor is limited, with SR8278 being the only identified synthetic antagonist. Moreover, no X-ray crystal structures are currently available that demonstrate the binding of REV-ERBα to antagonist ligands. This lack of structural information impedes the development of targeted therapeutics. To address this issue, we employed Gaussian accelerated molecular dynamics (GaMD) simulations to investigate the binding pathway of SR8278 to REV-ERBα. For comparison, we also used GaMD to observe the ligand binding process of STL1267, for which an X-ray structure is available. GaMD simulations successfully captured the binding of both ligands to the receptor's orthosteric site and predicted the ligand binding pathway and important amino acid residues involved in the antagonist SR8278 binding. This study highlights the effectiveness of GaMD in investigating protein-ligand interactions, particularly in the context of drug recognition for nuclear hormone receptors.

The switch states of the GDP-bound HRAS affected by point mutations: a study from Gaussian accelerated molecular dynamics simulations and free energy landscapes

Point mutations play a vital role in the conformational transformation of HRAS. In this work, Gaussian accelerated molecular dynamics (GaMD) simulations followed by constructions of free energy landscapes (FELs) were adopted to explore the effect of mutations D33K, A59T and L120A on conformation states of the GDP-bound HRAS. The results from the post-processing analyses on GaMD trajectories suggest that mutations alter the flexibility and motion modes of the switch domains from HRAS. The analyses from FELs show that mutations induce more disordered states of the switch domains and affect interactions of GDP with HRAS, implying that mutations yield a vital effect on the binding of HRAS to effectors. The GDP-residue interaction network revealed by our current work indicates that salt bridges and hydrogen bonding interactions (HBIs) play key roles in the binding of GDP to HRAS. Furthermore, instability in the interactions of magnesium ions and GDP with the switch SI leads to the extreme disorder of the switch domains. This study is expected to provide the energetic basis and molecular mechanism for further understanding the function of HRAS.Communicated by Ramaswamy H. Sarma.

Flanking Domains Modulate α-Synuclein Monomer Structure: A Molecular Dynamics Domain Deletion Study

Aggregates of misfolded α-synuclein proteins (asyn) are key markers of Parkinson's disease. Asyn proteins have three domains: an N-terminal domain, a hydrophobic NAC core implicated in aggregation, and a proline-rich C-terminal domain. Proteins with truncated C-terminal domains are known to be prone to aggregation and suggest that studying domain-domain interactions in asyn monomers could help elucidate the role of the flanking domains in modulating protein structure. To this end, we used Gaussian accelerated molecular dynamics (GAMD) to simulate wild-type (WT), N-terminal truncated (DN), C-terminal truncated (ΔC), and isolated NAC domain variants (isoNAC). Using clustering and contact analysis, we found that N- and C-terminal domains interact via electrostatic interactions, while the NAC and N-terminal domains interact through hydrophobic contacts. Our work also suggests that the C-terminal domain does not interact directly with the NAC domain but instead interacts with the N-terminal domain. Removal of the N-terminal domain led to increased contacts between NAC and C-terminal domains and the formation of interdomain β-sheets. Removal of either flanking domain also resulted in increased compactness of every domain. We also found that the contacts between flanking domains results in an electrostatic potential (ESP) that could possibly lead to favorable interactions with anionic lipid membranes. Removal of the C-terminal domain disrupts the ESP in a way that is likely to over-stabilize protein-membrane interactions. All of this suggests that one of the roles of the flanking domains may be to modulate the protein structure in a way that helps maintain elongation, hide hydrophobic residue from the solvent, and maintain an ESP that aids favorable interactions with the membrane.

Mutations influence the conformational dynamics of the GDP/KRAS complex

Mutations near allosteric sites can have a significant impact on the function of KRAS. Three specific mutations, K104Q, G12D/K104Q, and G12D/G75A, which are located near allosteric positions, were selected to investigate the molecular mechanisms behind mutation-induced influences on the activity of KRAS. Gaussian accelerated molecular dynamics (GaMD) simulations followed by the principal component analysis (PCA) were performed to improve the sampling of conformational states. The results revealed that these mutations significantly alter the structural flexibility, correlated motions, and dynamic behavior of the switch regions that are essential for KRAS binding to effectors or regulators. Furthermore, the mutations have a significant impact on the hydrogen bonding interactions between GDP and the switch regions, as well as on the electrostatic interactions of magnesium ions (Mg2+) with these regions. Our results verified that these mutations strongly influence the binding of KRAS to its effectors or regulators and allosterically regulate the activity. We believe that this work can provide valuable theoretical insights into a deeper understanding of KRAS function.Communicated by Ramaswamy H. Sarma.

Multiscale Study on the Intramolecular C-S Bond Formation Catalyzed by P450 Monooxygenase CxnD Involved in the Biosynthesis of Chuangxinmycin: The Critical Roles of Noncrystal Water Molecule and Conformational Change

Cytochrome P450 monooxygenase CxnD catalyzes intramolecular C-S bond formation in the biosynthesis of chuangxinmycin, which is representative of the synthesis of sulfur-containing natural heterocyclic compounds. The intramolecular cyclization usually requires the activation of two reaction sites and a large conformational change; thus, illuminating its detailed reaction mechanism remains challengeable. Here, the reaction pathway of CxnD-catalyzed C-S bond formation was clarified by a series of calculations, including Gaussian accelerated molecular dynamics simulations and quantum mechanical-molecular mechanical calculations. Our results revealed that the C-S formation follows a diradical coupling mechanism. CxnD first employs Cpd I to abstract the hydrogen atom from the imino group of the indole ring, and then, the resulted Cpd II further extracts another hydrogen atom from the thiol group of the side chain to afford a diradical intermediate, in which a noncrystal water molecule entering into the active site after the formation of Cpd I was proved to play an indispensable role. Moreover, the diradical intermediate cannot directly perform the coupling reaction. It should first undergo a series of conformational changes leading to the proximity of two reaction sites. It is the flexibility of the active site of the enzyme and the side chain of the substrate that makes the diradical coupling to be successful.

Mechanistic insights into ligand dissociation from the SARS-CoV-2 spike glycoprotein

The COVID-19 pandemic, driven by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has spurred an urgent need for effective therapeutic interventions. The spike glycoprotein of the SARS-CoV-2 is crucial for infiltrating host cells, rendering it a key candidate for drug development. By interacting with the human angiotensin-converting enzyme 2 (ACE2) receptor, the spike initiates the infection of SARS-CoV-2. Linoleate is known to bind the spike glycoprotein, subsequently reducing its interaction with ACE2. However, the detailed mechanisms underlying the protein-ligand interaction remain unclear. In this study, we characterized the pathways of ligand dissociation and the conformational changes associated with the spike glycoprotein by using ligand Gaussian accelerated molecular dynamics (LiGaMD). Our simulations resulted in eight complete ligand dissociation trajectories, unveiling two distinct ligand unbinding pathways. The preference between these two pathways depends on the gate distance between two α-helices in the receptor binding domain (RBD) and the position of the N-linked glycan at N343. Our study also highlights the essential contributions of K417, N121 glycan, and N165 glycan in ligand unbinding, which are equally crucial in enhancing spike-ACE2 binding. We suggest that the presence of the ligand influences the motions of these residues and glycans, consequently reducing accessibility for spike-ACE2 binding. These findings enhance our understanding of ligand dissociation from the spike glycoprotein and offer significant implications for drug design strategies in the battle against COVID-19.

Familial Alzheimer mutations stabilize synaptotoxic γ-secretase-substrate complexes

Mutations that cause familial Alzheimer's disease (FAD) are found in amyloid precursor protein (APP) and presenilin, the catalytic component of γ-secretase, that together produce amyloid β-peptide (Aβ). Nevertheless, whether Aβ is the primary disease driver remains controversial. We report here that FAD mutations disrupt initial proteolytic events in the multistep processing of APP substrate C99 by γ-secretase. Cryoelectron microscopy reveals that a substrate mimetic traps γ-secretase during the transition state, and this structure aligns with activated enzyme-substrate complex captured by molecular dynamics simulations. In silico simulations and in cellulo fluorescence microscopy support stabilization of enzyme-substrate complexes by FAD mutations. Neuronal expression of C99 and/or presenilin-1 in Caenorhabditis elegans leads to synaptic loss only with FAD-mutant transgenes. Designed mutations that stabilize the enzyme-substrate complex and block Aβ production likewise led to synaptic loss. Collectively, these findings implicate the stalled process-not the products-of γ-secretase cleavage of substrates in FAD pathogenesis.

Deciphering the molecular choreography of Janus kinase 2 inhibition via Gaussian accelerated molecular dynamics simulations: a dynamic odyssey

The Janus kinases (JAK) are crucial targets in drug development for several diseases. However, accounting for the impact of possible structural rearrangements on the binding of different kinase inhibitors is complicated by the extensive conformational variability of their catalytic kinase domain (KD). The dynamic KD contains mainly four prominent mobile structural motifs: the phosphate-binding loop (P-loop), the αC-helix within the N-lobe, the Asp-Phe-Gly (DFG) motif, and the activation loop (A-loop) within the C-lobe. These distinct structural orientations imply a complex signal transmission path for regulating the A-loop's flexibility and conformational preference for optimal JAK function. Nevertheless, the precise dynamical features of the JAK induced by different types of inhibitors still remain elusive. We performed comparative, microsecond-long, Gaussian accelerated molecular dynamics simulations in triplicate of three phosphorylated JAK2 systems: the KD alone, type-I ATP-competitive inhibitor (CI) bound KD in the catalytically active DFG-in conformation, and the type-II inhibitor (AI) bound KD in the catalytically inactive DFG-out conformation. Our results indicate significant conformational variations observed in the A-loop and αC helix motions upon inhibitor binding. Our studies also reveal that the DFG-out inactive conformation is characterized by the closed A-loop rearrangement, open catalytic cleft of N and C-lobe, the outward movement of the αC helix, and open P-loop states. Moreover, the outward positioning of the αC helix impacts the hallmark salt bridge formation between Lys882 and Glu898 in an inactive conformation. Finally, we compared their ligand binding poses and free energy by the MM/PBSA approach. The free energy calculations suggested that the AI's binding affinity is higher than CI against JAK2 due to an increased favorable contribution from the total non-polar interactions and the involvement of the αC helix. Overall, our study provides the structural and energetic insights crucial for developing more promising type I/II JAK2 inhibitors for treating JAK-related diseases.

Conformational States of the GDP- and GTP-Bound HRAS Affected by A59E and K117R: An Exploration from Gaussian Accelerated Molecular Dynamics

The HRAS protein is considered a critical target for drug development in cancers. It is vital for effective drug development to understand the effects of mutations on the binding of GTP and GDP to HRAS. We conducted Gaussian accelerated molecular dynamics (GaMD) simulations and free energy landscape (FEL) calculations to investigate the impacts of two mutations (A59E and K117R) on GTP and GDP binding and the conformational states of the switch domain. Our findings demonstrate that these mutations not only modify the flexibility of the switch domains, but also affect the correlated motions of these domains. Furthermore, the mutations significantly disrupt the dynamic behavior of the switch domains, leading to a conformational change in HRAS. Additionally, these mutations significantly impact the switch domain's interactions, including their hydrogen bonding with ligands and electrostatic interactions with magnesium ions. Since the switch domains are crucial for the binding of HRAS to effectors, any alterations in their interactions or conformational states will undoubtedly disrupt the activity of HRAS. This research provides valuable information for the design of drugs targeting HRAS.

Biased Signaling in Mutated Variants of β2-Adrenergic Receptor: Insights from Molecular Dynamics Simulations

The molecular basis of receptor bias in G protein-coupled receptors (GPCRs) caused by mutations that preferentially activate specific intracellular transducers over others remains poorly understood. Two experimentally identified biased variants of β2-adrenergic receptors (β2AR), a prototypical GPCR, are a triple mutant (T68F, Y132A, and Y219A) and a single mutant (Y219A); the former bias the receptor toward the β-arrestin pathway by disfavoring G protein engagement, while the latter induces G protein signaling explicitly due to selection against GPCR kinases (GRKs) that phosphorylate the receptor as a prerequisite of β-arrestin binding. Though rigorous characterizations have revealed functional implications of these mutations, the atomistic origin of the observed transducer selectivity is not clear. In this study, we investigated the allosteric mechanism of receptor bias in β2AR using microseconds of all-atom Gaussian accelerated molecular dynamics (GaMD) simulations. Our observations reveal distinct rearrangements in transmembrane helices, intracellular loop 3, and critical residues R1313.50 and Y3267.53 in the conserved motifs D(E)RY and NPxxY for the mutant receptors, leading to their specific transducer interactions. Moreover, partial dissociation of G protein from the receptor core is observed in the simulations of the triple mutant in contrast to the single mutant and wild-type receptor. The reorganization of allosteric communications from the extracellular agonist BI-167107 to the intracellular receptor-transducer interfaces drives the conformational rearrangements responsible for receptor bias in the single and triple mutants. The molecular insights into receptor bias of β2AR presented here could improve the understanding of biased signaling in GPCRs, potentially opening new avenues for designing novel therapeutics with fewer side-effects and superior efficacy.

Role of ATP Hydrolysis and Product Release in the Translocation Mechanism of SARS-CoV-2 NSP13

In response to the emergence of COVID-19, caused by SARS-CoV-2, there has been a growing interest in understanding the functional mechanisms of the viral proteins to aid in the development of new therapeutics. Nonstructural protein 13 (nsp13) helicase is an attractive target for antivirals because it is essential for viral replication and has a low mutation rate, yet the structural mechanisms by which this enzyme binds and hydrolyzes ATP to cause unidirectional RNA translocation remain elusive. Using Gaussian accelerated molecular dynamics (GaMD), we generated comprehensive conformational ensembles of all substrate states along the ATP-dependent cycle. Shape-GMM clustering of the protein yields four protein conformations that describe an opening and closing of both the ATP pocket and the RNA cleft that is achieved through a combination of conformational selection and induction along the ATP hydrolysis cycle. Furthermore, three protein-RNA conformations are observed that implicate motifs Ia, IV, and V as playing a pivotal role in an ATP-dependent inchworm translocation mechanism. Finally, based on a linear discriminant analysis of protein conformations, we identify L405 as a pivotal residue for the opening and closing mechanism and propose a L405D mutation as a way to disrupt translocation. This research enhances our understanding of nsp13's role in viral replication and could contribute to the development of antiviral strategies.

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.

Characterization of Posttranslationally Modified PHF-1 Tau Peptides Using Gaussian Accelerated Molecular Dynamics Simulation

The microtubule-associated protein, Tau, is an intrinsically disordered protein that plays a crucial role in neurodegenerative diseases like Alzheimer's disease. The posttranslational modifications across the Tau protein domains are involved in regulating Tau protein's function and disease onset. Of the various posttranslational modifications at Ser, Thr, and Tyr sites, O-GlcNAcylation and phosphorylation are the most critical ones, playing a vital role in Tau aggregation and tauopathies. To understand the function, it is essential to characterize the structural changes associated with Tau modification. Previous experimental studies have focused on high-resolution nuclear magnetic resonance techniques to structurally characterize the effect of phosphorylation, O-GlcNAcylation, and combination of both PTMs on Tau conformation in small peptides centered on the PHF-1 epitope from amino acid 392 to 411. The structural characterization using atomistic molecular dynamics simulation of such disordered peptides requires long simulation time, proper sampling method, and utilization of appropriate force fields for accurate determination of conformational ensembles, resembling the experimental data. This chapter details the protocol for the structural characterization of modified Tau peptides using the CHARMM36m force field and enhanced sampling methods like Gaussian accelerated molecular dynamics (GaMD) simulation. We have focused on a detailed explanation of the GaMD method and analyses of molecular dynamics trajectories to explain the relationship between two modifications, phospho- and glyco-, at C-terminus of Tau protein and its stable conformation over the longer simulation timeframes. The analyses involve energetics reweighting, clustering of simulation trajectories, and characterization of secondary structure using circular dichroism data from the simulation. The reader can utilize this protocol to investigate the structures of complex proteins, especially the disordered ones.

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.

When MELD Meets GaMD: Accelerating Biomolecular Landscape Exploration

We introduce Gaussian accelerated MELD (GaMELD) as a new method for exploring the energy landscape of biomolecules. GaMELD combines the strengths of Gaussian accelerated molecular dynamics (GaMD) and modeling employing limited data (MELD) to navigate complex energy landscapes. MELD uses replica-exchange molecular simulations to integrate limited and uncertain data into simulations via Bayesian inference. MELD has been successfully applied to problems of structure prediction like protein folding and complex structure prediction. However, the computational cost for MELD simulations has limited its broader applicability. The synergy of GaMD and MELD surmounts this limitation efficiently sampling the energy landscape at a lower computational cost (reducing the computational cost by a factor of 2 to six). Effectively, GaMD is used to shift energy distributions along replicas to increase the overlap in energy distributions across replicas, facilitating a random walk in replica space. We tested GaMELD on a benchmark set of 12 small proteins that have been previously studied through MELD and conventional MD. GaMELD consistently achieves accurate predictions with fewer replicas. By increasing the efficacy of replica exchange, GaMELD effectively accelerates convergence in the conformational space, enabling improved sampling across a diverse set of systems.

Molecular Dynamics Activation of γ-Secretase for Cleavage of the Notch1 Substrate

γ-Secretase is an intramembrane aspartyl protease complex that cleaves the transmembrane domain of over 150 peptide substrates, including amyloid precursor protein (APP) and the Notch family of receptors, via two conserved aspartates D257 and D385 in the presenilin-1 (PS1) catalytic subunit. However, while the activation of γ-secretase for cleavage of APP has been widely studied, the cleavage of Notch by γ-secretase remains poorly explored. Here, we combined Gaussian accelerated molecular dynamics (GaMD) simulations and mass spectrometry (MS) analysis of proteolytic products to present the first dynamic models for cleavage of Notch by γ-secretase. MS showed that γ-secretase cleaved the WT Notch at Notch residue G34, while cleavage of the L36F mutant Notch occurred at Notch residue C33. Initially, we prepared our simulation systems starting from the cryoEM structure of Notch-bound γ-secretase (PDB: 6IDF) and failed to capture the proper cleavages of WT and L36F Notch by γ-secretase. We then discovered an incorrect registry of the Notch substrate in the PS1 active site through alignment of the experimental structure of Notch-bound (PDB: 6IDF) and APP-bound γ-secretase (PDB: 6IYC). Every residue of the APP substrate was systematically mutated to the corresponding Notch residue to prepare a resolved model of Notch-bound γ-secretase complexes. GaMD simulations of the resolved model successfully captured γ-secretase activation for proper cleavages of both WT and L36F mutant Notch. Our findings presented here provided mechanistic insights into the structural dynamics and enzyme-substrate interactions required for γ-secretase activation for cleavage of Notch and other substrates.

An in silico investigation of Kv2.1 potassium channel: Model building and inhibitors binding sites analysis

Kv2.1 is widely expressed in brain, and inhibiting Kv2.1 is a potential strategy to prevent cell death and achieve neuroprotection in ischemic stroke. Herein, an in silico model of Kv2.1 tetramer structure was constructed by employing the AlphaFold-Multimer deep learning method to facilitate the rational discovery of Kv2.1 inhibitors. GaMD was utilized to create an ion transporting trajectory, which was analyzed with HMM to generate multiple representative receptor conformations. The binding site of RY785 and RY796(S) under the P-loop was defined with Fpocket program together with the competitive binding electrophysiology assay. The docking poses of the two inhibitors were predicted with the aid of the semi-empirical quantum mechanical calculation, and the IGMH results suggested that Met375, Thr376, and Thr377 of the P-helix and Ile405 of the S6 segment made significant contributions to the binding affinity. These results provided insights for rational molecular design to develop novel Kv2.1 inhibitors.

Predicting absolute aqueous solubility by applying a machine learning model for an artificially liquid-state as proxy for the solid-state

In this study, we use machine learning algorithms with QM-derived COSMO-RS descriptors, along with Morgan fingerprints, to predict the absolute solubility of drug-like compounds. The QM-derived descriptors account for the molecular properties of the solute, i.e., the solute-solute interactions in an artificial-liquid-state (super-cooled liquid), and the solute-solvent interactions in solution. We employ two main approaches to predict solubility: (i) a hypothetical pathway that involves melting the solute at room temperature T = T¯ ([Formula: see text]) and mixing the artificially liquid solute into the solvent ([Formula: see text]). In this approach [Formula: see text] is predicted using machine learning models, and the [Formula: see text] is obtained from COSMO-RS calculations; (ii) direct solubility prediction using machine learning algorithms. The models were trained on a large number of Bayer in-house compounds for which water solubility data is available at physiological pH of 6.5 and ambient temperature. We also evaluated our models using external datasets from a solubility challenge. Our models present great improvements compared to the absolute solubility prediction with the QSAR model for the artificial liquid state as implemented in the COSMOtherm software, for both in-house and external datasets. We are furthermore able to demonstrate the superiority of QM-derived descriptors compared to cheminformatics descriptors. We finally present low-cost alternative models using fragment-based COSMOquick calculations with only marginal reduction in the quality of predicted solubility.

Modeling the Effect on a Novel Fungal Peptaibol Placed in an All-Atom Bacterial Membrane Mimicking System via Accelerated Molecular Dynamics Simulations

We previously reported on a novel peptaibol, named Tripleurin XIIc (TPN), an 18-residue long sequence produced by the fungus Trichoderma pleuroti. We elucidated its 3D structure via classical and accelerated molecular dynamics simulation (aMD) methods and reported the folding dynamics of TPN in water and chloroform solvents. Peptaibols, in general, are insoluble in water, as they are amphipathic and may prefer hydrophobic environments like transmembrane regions. In this study, we attempted to use aMD simulations to model an all-atom bacterial membrane system while placing a TPN molecule in its vicinity. The results highlighted that TPN was able to introduce some disorder into the membrane and caused lipid clustering. It could also enter the transmembrane region from the water-bilayer interface. The structural dynamics of TPN in the transmembrane region revealed a single energetically stable conformation similar to the one obtained from water and chloroform solvent simulations reported by us previously. However, this linear structure was found to be at the local energy minimum (stable) in water but at a metastable intermediate state (higher energy) in chloroform. Therefore, it could be said that the water solvent can be successfully used for folding simulations of peptaibols.

Addressing Challenges of Macrocyclic Conformational Sampling in Polar and Apolar Solvents: Lessons for Chameleonicity

We evaluated a workflow to reliably sample the conformational space of a set of 47 peptidic macrocycles. Starting from SMILES strings, we use accelerated molecular dynamics simulations to overcome high energy barriers, in particular, the cis-trans isomerization of peptide bonds. We find that our approach performs very well in polar solvents like water and dimethyl sulfoxide. Interestingly, the protonation state of a secondary amine in the ring only slightly influences the conformational ensembles of our test systems. For several of the macrocycles, determining the conformational distribution in chloroform turns out to be considerably more challenging. Especially, the choice of partial charges crucially influences the ensembles in chloroform. We address these challenges by modifying initial structures and the choice of partial charges. Our results suggest that special care has to be taken to understand the configurational distribution in apolar solvents, which is a key step toward a reliable prediction of membrane permeation of macrocycles and their chameleonic properties.

Deep Learning Dynamic Allostery of G-Protein-Coupled Receptors

G-protein-coupled receptors (GPCRs) make up the largest superfamily of human membrane proteins and represent primary targets of ∼1/3 of currently marketed drugs. Allosteric modulators have emerged as more selective drug candidates compared with orthosteric agonists and antagonists. However, many X-ray and cryo-EM structures of GPCRs resolved so far exhibit negligible differences upon the binding of positive and negative allosteric modulators (PAMs and NAMs). The mechanism of dynamic allosteric modulation in GPCRs remains unclear. In this work, we have systematically mapped dynamic changes in free energy landscapes of GPCRs upon binding of allosteric modulators using the Gaussian accelerated molecular dynamics (GaMD), deep learning (DL), and free energy prOfiling Workflow (GLOW). GaMD simulations were performed for a total of 66 μs on 44 GPCR systems in the presence and absence of the modulator. DL and free energy calculations revealed significantly reduced dynamic fluctuations and conformational space of GPCRs upon modulator binding. While the modulator-free GPCRs often sampled multiple low-energy conformational states, the NAMs and PAMs confined the inactive and active agonist-G-protein-bound GPCRs, respectively, to mostly only one specific conformation for signaling. Such cooperative effects were significantly reduced for binding of the selective modulators to "non-cognate" receptor subtypes. Therefore, GPCR allostery exhibits a dynamic "conformational selection" mechanism. In the absence of available modulator-bound structures as for most current GPCRs, it is critical to use a structural ensemble of representative GPCR conformations rather than a single structure for compound docking ("ensemble docking"), which will potentially improve structure-based design of novel allosteric drugs of GPCRs.

Molecular analysis of the type III interferon complex and its applications in protein engineering

Type III interferons (IFNλs) are cytokines with critical roles in the immune system and are attractive therapeutic candidates due to their tissue-specific activity. Despite entering several clinical trials, results have demonstrated limited efficacy and potency, partially attributed to low-affinity protein-protein interactions (PPIs) responsible for receptor complex formation. Subsequently, structural studies of the native IFNλ signaling complexes remain inaccessible. While protein engineering can overcome affinity limitations, tools to investigate low-affinity systems like these remain limited. To provide insights into previous efforts to strengthen the PPIs within this complex, we perform a molecular analysis of the extracellular ternary complexes of IFNλ3 using both computational and experimental approaches. We first use molecular simulations and modeling to quantify differences in PPIs and residue strain fluctuations, generate detailed free energy landscapes, and reveal structural differences between an engineered, high-affinity complex, and a model of the wild-type, low-affinity complex. This analysis illuminates distinct behaviors of these ligands, yielding mechanistic insights into IFNλ complex formation. We then apply these computational techniques in protein engineering and design by utilizing simulation data to identify hotspots of interaction to rationally engineer the native cytokine-receptor complex for increased stability. These simulations are then validated by experimental techniques, showing that a single mutation at a computationally predicted site of interaction between the two receptors increases PPIs and improves complex formation for all IFNλs. This study highlights the power of molecular dynamics simulations for protein engineering and design as applied to the IFNλ family but also presents a potential tool for analysis and engineering of other systems with low-affinity PPIs.

The role of glycerol-water mixtures in the stability of FKBP12-rapalog-FRB complexes

The thermodynamic and biophysical implications of the introduction of a co-solvent during protein-ligand binding remain elusive. Using ternary complexes of 12-kDa FK506 binding protein (FKBP12), FKBP-rapamycin binding (FRB) domain of the mammalian/mechanistic target of rapamycin (mTOR) kinase, and rapamycin analogs (rapalogs) in glycerol-water mixtures, the influence of solvent composition on ligand binding dynamics was explored. The pharmaceutical potential of rapalogs and the utility of glycerol as a co-solvent in drug delivery applications were critical in deciding the system to be studied. Consolidation of existing studies on rapamycin modification was first performed to strategically design a new rapalog called T1. The results from 100-ns dual-boost Gaussian accelerated molecular dynamics simulations showed that protein stability was induced in the presence of glycerol. Reweighting of the trajectories revealed that the glycerol-rich solvent system lowers the energy barrier in the conformational space of the protein while also preserving native contacts between the ligand and the residues in the binding site. Calculated binding free energies using MM/GBSA also showed that electrostatic energy and polar contribution of solvation energy are heavily influenced by the changes in solvation. Glycerol molecules are preferentially excluded through electrostatic interactions from the solvation shell which induce complex stability as seen in existing experiments. Hence, using glycerol as a co-solvent in rapamycin delivery has a significant role in maintaining stability. In addition, compound T1 is a potential mTORC1-selective inhibitor with strong affinity for the FKBP12-FRB complex. This study aims to provide insights on the design of new rapalogs, and the applicability of glycerol as co-solvent for FKBP12-rapalog-FRB complexes.

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.

Dynamic Insights into the Self-Activation Pathway and Allosteric Regulation of the Orphan G-Protein-Coupled Receptor GPR52

Within over 800 members of G-protein-coupled receptors, there are numerous orphan receptors whose endogenous ligands are largely unknown, providing many opportunities for novel drug discovery. However, the lack of an in-depth understanding of the intrinsic working mechanism for orphan receptors severely limits the related rational drug design. The G-protein-coupled receptor 52 (GPR52) is a unique orphan receptor that constitutively increases cellular 5'-cyclic adenosine monophosphate (cAMP) levels without binding any exogenous agonists and has been identified as a promising therapeutic target for central nervous system disorders. Although recent structural biology studies have provided snapshots of both active and inactive states of GPR52, the mechanism of the conformational transition between these states remains unclear. Here, an acceptable self-activation pathway for GPR52 was proposed through 6 μs Gaussian accelerated molecular dynamics (GaMD) simulations, in which the receptor spontaneously transitions from the active state to that matching the inactive crystal structure. According to the three intermediate states of the receptor obtained by constructing a reweighted potential of mean force, how the allosteric regulation occurs between the extracellular orthosteric binding pocket and the intracellular G-protein-binding site is revealed. Combined with the independent gradient model, several important microswitch residues and the allosteric communication pathway that directly links the two regions are both identified. Transfer entropy calculations not only reveal the complex allosteric signaling within GPR52 but also confirm the unique role of ECL2 in allosteric regulation, which is mutually validated with the results of GaMD simulations. Overall, this work elucidates the allosteric mechanism of GPR52 at the atomic level, providing the most detailed information to date on the self-activation of the orphan receptor.

Multiscale computational modeling of the effects of 2'-deoxy-ATP on cardiac muscle calcium handling

2'-Deoxy-ATP (dATP), a naturally occurring near analog of ATP, is a well-documented myosin activator that has been shown to increase contractile force, improve pump function, and enhance lusitropy in the heart. Calcium transients in cardiomyocytes with elevated levels of dATP show faster calcium decay compared with cardiomyocytes with basal levels of dATP, but the mechanisms behind this are unknown. Here, we design and utilize a multiscale computational modeling framework to test the hypothesis that dATP acts on the sarcoendoplasmic reticulum calcium-ATPase (SERCA) pump to accelerate calcium re-uptake into the sarcoplasmic reticulum during cardiac relaxation. Gaussian accelerated molecular dynamics simulations of human cardiac SERCA2A in the E1 apo, ATP-bound and dATP-bound states showed that dATP forms more stable contacts in the nucleotide binding pocket of SERCA and leads to increased closure of cytosolic domains. These structural changes ultimately lead to changes in calcium binding, which we assessed using Brownian dynamics simulations. We found that dATP increases calcium association rate constants to SERCA and that dATP binds to apo SERCA more rapidly than ATP. Using a compartmental ordinary differential equation model of human cardiomyocyte excitation-contraction coupling, we found that these increased association rate constants contributed to the accelerated rates of calcium transient decay observed experimentally. This study provides clear mechanistic evidence of enhancements in cardiac SERCA2A pump function due to interactions with dATP.

Insights into dynamic changes in ADC-7 and P99 cephalosporinases using small angle x-ray scattering (SAXS)

To counter the emergence of β-lactamase (BL) mediated resistance, design of new β-lactamase inhibitors (BLIs) is critical. Many high-resolution crystallographic structures of BL complexed with BLIs are available. However, their impact on BLI design is struggling to keep pace with novel and emerging variants. Small angle x-ray scattering (SAXS) in combination with molecular modeling is a useful tool to determine dynamic structures of macromolecules in solution. An important application of SAXS is to determine the conformational changes that occur when BLI bind to BL. To probe if conformational dynamics occur in class C cephalosporinases, we studied SAXS profiles of two clinically relevant class C β-lactamases, Acinetobacter baumannii ADC-7 and Enterobacter cloacae P99 in apo format complexed with BLIs. Importantly, SAXS data analysis demonstrated that in solution, these representative class C enzymes remain monomeric and did not show the associated assemblies that were seen in various crystal structures. SAXS data acquired for ADC-7 and P99, in apo and inhibitor bound states, clearly showed that these enzymes undergo detectable conformational changes, and these class C β-lactamases also close upon binding inhibitors as does BlaC. Further analysis revealed that addition of inhibitor led to the compacting of a range of residues around the active site, indicating that the conformational changes that both P99 and ADC-7 undergo are central to inhibitor recognition and efficacy. Our findings support the importance of exploring conformational changes using SAXS analysis in the design of future BLIs.Communicated by Ramaswamy H. Sarma.

Enhanced Configurational Sampling Approaches to Alchemical Ligand Binding Free Energy Simulations: Current Status and Challenges

Ligand binding free energy simulations (LB-FES) have been routine tasks in modern drug discovery campaign. A long-standing challenge for LB-FES is the difficulty in adequately sampling nontrivial environmental reorganizations in response to ligand binding. Therefore, various enhanced configurational sampling (ECS) approaches were devised to speed up fluctuations of relevant slow degrees of freedom (SDOF) and ensure simulation convergence. However, in contrast to the achievements in parametrization, software performance, and workflow automation, efficient ECS methodology suitable for high throughput screening remains in an early stage of development. Here, a review of ECS developments with LB-FES is presented, revisiting current approaches and underlining the major technical pitfalls and challenges. This Perspective focuses on alchemical LB-FES on account of their predominant role in high throughput drug screening as well as the established partnership with ECS. The critical aspects of designing ECS approaches, from both theoretical and applied perspectives, are described. This work is intended to provide a contemporary review of the scientific, technical, and practical issues associated with the accelerating convergence of alchemical LB-FES.

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.

Structure-Based Discovery of a Novel Allosteric Inhibitor against Human Dopamine Transporter

Human dopamine transporter (hDAT) regulates the reuptake of extracellular dopamine (DA) and is an essential therapeutic target for central nervous system (CNS) diseases. The allosteric modulation of hDAT has been identified for decades. However, the molecular mechanism underlying the transportation is still elusive, which hinders the rational design of allosteric modulators against hDAT. Here, a systematic structure-based method was performed to explore allosteric sites on hDAT in inward-open (IO) conformation and to screen compounds with allosteric affinity. First, the model of the hDAT structure was constructed based on the recently reported Cryo-EM structure of the human serotonin transporter (hSERT) and Gaussian-accelerated molecular dynamics (GaMD) simulation was further utilized for the identification of intermediate energetic stable states of the transporter. Then, with the potential druggable allosteric site on hDAT in IO conformation, virtual screening of seven enamine chemical libraries (∼440,000 compounds) was processed, resulting in 10 compounds being purchased for in vitro assay and with Z1078601926 discovered to allosterically inhibit hDAT (IC₅₀ = 0.527 [0.284; 0.988] μM) when nomifensine was introduced as an orthosteric ligand. Finally, the synergistic effect underlying the allosteric inhibition of hDAT by Z1078601926 and nomifensine was explored using additional GaMD simulation and postbinding free energy analysis. The hit compound discovered in this work not only provides a good starting point for lead optimization but also demonstrates the usability of the method for the structure-based discovery of novel allosteric modulators of other therapeutic targets.

Extended Ensemble Molecular Dynamics Study of Ammonia-Cellulose I Complex Crystal Models: Free-Energy Landscape and Atomistic Pictures of Ammonia Diffusion in the Crystalline Phase

Here, we report extended ensemble molecular dynamics simulations of ammonia-cellulose I complex crystal models to evaluate the diffusion behavior of the guest ammonia molecules and the potential of mean force (PMF), namely, the free energy change along the chosen reaction coordinate, for migration of an ammonia molecule in the crystal models. Accelerated molecular dynamics simulations confirmed that ammonia molecules almost exclusively diffused through the hydrophilic channel even when the crystal framework was retained. Adaptive steered molecular dynamics simulations detected distinct PMF peaks with heights of approximately 7 kcal/mol as the ammonia molecule passed through the cellulose-chain layers. Introducing hybrid quantum mechanical and molecular mechanics theory to the adaptive steered molecular dynamics simulation effectively lowered the heights of the PMF peaks to approximately 5 kcal/mol, accompanied by a slight decrease in the baseline. Removal of the ammonia molecules in the neighboring channels resulted in a continuous increase in the baseline for the migration of an ammonia molecule in the hydrophilic channel. When the halves of the crystal model were separated to widen the hydrophilic channel to 0.2 nm, the PMF profiles exhibited an unexpected increase. This resulted from water structuring in the expanded hydrophilic channel, which disappeared with further expansion of the hydrophilic channel to 0.3 nm.

A New Approach for Estimating the Free Energy Differences among Multiple Thermodynamic States in Statistical Simulations

In this letter, a new approach to compute free energy differences (FEDs) between multiple thermodynamics states is introduced. The method directly uses energy probability densities, which can be extracted with high accuracy from equilibrium simulations to obtain FEDs. Methods in current use, such as Bennett acceptance ratio (BAR), its multistate generalization (MBAR), or the weighted histogram analysis method (WHAM), require iterative solution of nonlinear equations which are known to be slowly convergent. The equations providing MBAR FEDs are identical to those derived earlier by Souaille and Roux in a method that has become known informally as "binless WHAM". In contrast, we obtain FEDs by solution of linear equations. For the classic two-state problem, the statistical error of our method, solving linear equations, is shown analytically to match that of BAR under common conditions.

Pharmacological hallmarks of allostery at the M4 muscarinic receptor elucidated through structure and dynamics

Allosteric modulation of G protein-coupled receptors (GPCRs) is a major paradigm in drug discovery. Despite decades of research, a molecular-level understanding of the general principles that govern the myriad pharmacological effects exerted by GPCR allosteric modulators remains limited. The M₄ muscarinic acetylcholine receptor (M₄ mAChR) is a validated and clinically relevant allosteric drug target for several major psychiatric and cognitive disorders. In this study, we rigorously quantified the affinity, efficacy, and magnitude of modulation of two different positive allosteric modulators, LY2033298 (LY298) and VU0467154 (VU154), combined with the endogenous agonist acetylcholine (ACh) or the high-affinity agonist iperoxo (Ipx), at the human M₄ mAChR. By determining the cryo-electron microscopy structures of the M₄ mAChR, bound to a cognate G_i1 protein and in complex with ACh, Ipx, LY298-Ipx, and VU154-Ipx, and applying molecular dynamics simulations, we determine key molecular mechanisms underlying allosteric pharmacology. In addition to delineating the contribution of spatially distinct binding sites on observed pharmacology, our findings also revealed a vital role for orthosteric and allosteric ligand-receptor-transducer complex stability, mediated by conformational dynamics between these sites, in the ultimate determination of affinity, efficacy, cooperativity, probe dependence, and species variability. There results provide a holistic framework for further GPCR mechanistic studies and can aid in the discovery and design of future allosteric drugs.

Accelerated Molecular Dynamics for Peptide Folding: Benchmarking Different Combinations of Force Fields and Explicit Solvent Models

Accelerated molecular dynamics (aMD) protocols were assessed on predicting the secondary structure of eight peptides, of which two are helical, three are β-hairpins, and three are disordered. Protocols consisted of combinations of three force fields (ff99SB, ff14SB, ff19SB) and two explicit solvation models (TIP3P and OPC), and were evaluated in two independent aMD simulations, one starting from an extended conformation, the other starting from a misfolded conformation. The results of these analyses indicate that all three combinations performed well on helical peptides. As for β-hairpins, ff19SB performed well with both solvation methods, with a slight preference for the TIP3P solvation model, even though performance was dependent on both peptide sequence and initial conformation. The ff19SB/OPC combination had the best performance on intrinsically disordered peptides. In general, ff14SB/TIP3P suffered the strongest helical bias.

Multi-state Model-Based Identification of Cryptic Allosteric Sites on Human Serotonin Transporter

Serotonin transporter (SERT) plays a fundamental role in taking the synaptic cleft serotonin back to the presynaptic neuron. The discovery of allosteric SERT modulators represents the next-generation medication for psychiatric disorders such as depression. Here, based on the cryo-EM structures of ibogaine in complex with SERT in distinct conformations, the multiple functional structures of the transporter bound to serotonin, including outward-open (OO_holo), outward-occluded (OC_holo), and inward-open (IO_holo and IO_holo'), were carefully characterized by induced-fit docking Gaussian-accelerated molecular dynamics (IFD-GaMD) simulation and the free-energy landscape analysis. Further MM/GBSA binding free energy, per-residue contribution, and molecular interaction fingerprint calculations revealed the interaction variations of serotonin with SERT in functional structures, which confirmed the allostery of SERT during serotonin reuptake. Moreover, five unique cryptic allosteric sites, which are druggable and capable of targeting by small molecules, were identified on the characterized multistate structures. These results provide structural and energetic information for the molecular mechanism of serotonin reuptake and will provide opportunities for the development of novel therapeutics based on the identified new allosteric sites on SERT.

Computational insights into ligand-induced G protein and β-arrestin signaling of the dopamine D1 receptor

The dopamine D1 receptor (D1R), is a class A G protein coupled-receptor (GPCR) which has been a promising drug target for psychiatric and neurological disorders such as Parkinson's disease (PD). Previous studies have suggested that therapeutic effects can be realized by targeting the β-arrestin signaling pathway of dopamine receptors, while overactivation of the G protein-dependent pathways leads to side effects, such as dyskinesias. Therefore, it is highly desirable to develop a D1R ligand that selectively regulates the β-arrestin pathway. Currently, most D1R agonists are signaling-balanced and stimulate both G protein and β-arrestin pathways, with a few reports of G protein biased ligands. However, identification and characterization of β-arrestin biased D1R agonists has been a challenge thus far. In this study, we implemented Gaussian accelerated molecular dynamics (GaMD) simulations to provide valuable computational insights into the possible underlying molecular mechanism of the different signaling properties of two catechol and two non-catechol D1R agonists that are either G protein biased or signaling-balanced. Dynamic network analysis further identified critical residues in the allosteric signaling network of D1R for each ligand at different conformational or binding states. Some of these residues are crucial for G protein or arrestin signals of GPCRs based on previous studies. Finally, we provided a molecular design strategy which can be utilized by medicinal chemists to develop potential β-arrestin biased D1R ligands. The proposed hypotheses are experimentally testable and can guide the development of safer and more effective medications for a variety of CNS disorders.

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/.

Impacts of Mutations in the P-Loop on Conformational Alterations of KRAS Investigated with Gaussian Accelerated Molecular Dynamics Simulations

G12 mutations heavily affect conformational transformation and activity of KRAS. In this study, Gaussian accelerated molecular dynamics (GaMD) simulations were performed on the GDP-bound wild-type (WT), G12A, G12D, and G12R KRAS to probe mutation-mediated impacts on conformational alterations of KRAS. The results indicate that three G12 mutations obviously affect the structural flexibility and internal dynamics of the switch domains. The analyses of the free energy landscapes (FELs) suggest that three G12 mutations induce more conformational states of KRAS and lead to more disordered switch domains. The principal component analysis shows that three G12 mutations change concerted motions and dynamics behavior of the switch domains. The switch domains mostly overlap with the binding region of KRAS to its effectors. Thus, the high disorder states and concerted motion changes of the switch domains induced by G12 mutations affect the activity of KRAS. The analysis of interaction network of GDP with KRAS signifies that the instability in the interactions of GDP and magnesium ion with the switch domain SW1 drives the high disordered state of the switch domains. This work is expected to provide theoretical aids for understanding the function of KRAS.

How the Conformational Movement of the Substrate Drives the Regioselective C-N Bond Formation in P450 TleB: Insights from Molecular Dynamics Simulations and Quantum Mechanical/Molecular Mechanical Calculations

P450 TleB catalyzes the oxidative cyclization of the dipeptide N-methylvalyl-tryptophanol into indolactam V through selective intramolecular C-H bond amination at the indole C4 position. Understanding its catalytic mechanism is instrumental for the engineering or design of P450-catalyzed C-H amination reactions. Using multiscale computational methods, we show that the reaction proceeds through a diradical pathway, involving a hydrogen atom transfer (HAT) from N1-H to Cpd I, a conformational transformation of the substrate radical species, and a second HAT from N13-H to Cpd II. Intriguingly, the conformational transformation is found to be the key to enabling efficient and selective C-N coupling between N13 and C4 in the subsequent diradical coupling reaction. The underlined conformational transformation is triggered by the first HAT, which proceeds with an energy-demanding indole ring flip and is followed by the facile approach of the N13-H group to Cpd II. Detailed analysis shows that the internal electric field (IEF) from the protein environment plays key roles in the transformation process, which not only provides the driving force but also stabilizes the flipped conformation of the indole radical. Our simulations provide a clear picture of how the P450 enzyme can smartly modulate the selective C-N coupling reaction. The present findings are in line with all available experimental data, highlighting the crucial role of substrate dynamics in controlling this highly valuable reaction.