Exploring ligand binding pathways on proteins using hypersound-accelerated molecular dynamics

Precision spatiotemporal analysis of large-scale compound-protein interactions through molecular dynamics simulation

Biological systems are composed of and regulated by intricate and diverse biomolecular interactions. Experimental and computational approaches have been developed to elucidate the mechanisms of these interactions; however, owing to cost, time, and accuracy issues, large-scale spatiotemporal analyses of molecular pairs remain challenging. Thus, the molecular recognition mechanisms underlying these diverse interactions remain unclear. We successfully simulated the large-scale molecular dynamics (MD) of 4,275 protein-compound pairs by combining a method to accelerate the MD simulations with the supercomputer Fugaku. Our spatiotemporal analysis of generated big MD data revealed universal features underlying molecular recognition and binding processes. This study expands our understanding of the concept of MD simulations from a technique to investigate the dynamic properties of individual protein-drug pairs to an approach to perform large-scale spatiotemporal analysis and compound screening. This study opens an avenue in biological research for subsequent drug discovery.

Pathways and mechanism of MRTX1133 binding to KRAS G12D elucidated by molecular dynamics simulations and Markov state models

KRAS G12D is the most common oncogenic mutation identified in several types of cancer. Therefore, design of inhibitors targeting KRAS G12D represents a promising strategy for anticancer therapy. MRTX1133 is a highly potent inhibitor (approximate experiment Kd ≈ 0.0002 nM) of KRAS G12D and is currently in Phase 1/2 study, however, pathways of the compound binding to KRAS G12D has remained unknown, and the mechanism underlying the complicated dynamic process are challenging to capture experimentally, which hinder the structure-based anti-cancer drug design. Here, using MRTX1133 as a probe, unbiased molecular dynamics (MD) was used to simulate the process of MRTX1133 spontaneously binding to KRAS G12D. In six of 42 independent MD simulation (a total of 99 μs), MRTX1133 was observed to successfully associate with KRAS G12D. The kinetically metastable states refer to the potential pathways of MRTX1133 binding to KRAS G12D were revealed by Markov state models (MSM) analysis. Additionally, 8 key residues that are essential for MRTX1133 recognition and tight binding at the preferred low energy states were identified by MM/GBSA analysis. In sum, this study provides a new perspective on understanding the pathways and mechanism of MRTX1133 binding to KRAS G12D.

Unveiling the Mechanisms of EGCG-p53 Interactions through Molecular Dynamics Simulations

Green tea consumption is associated with protective and preventive effects against various types of cancer. These effects are attributed to polyphenols, particularly epigallocatechin-3-gallate (EGCG). EGCG acts by directly inhibiting tumor suppressor protein p53. The binding mechanism by which EGCG inhibits p53 activity is associated with residues Trp23-Lys24 and Pro47-Thr55 within the p53 N-terminal domain (NTD). However, the structural and thermodynamic aspects of the interaction between EGCG and p53 are poorly understood. Therefore, based on crystallographic data, we combine docking, molecular dynamics (MD) simulations, and molecular mechanics generalized Born surface area approaches to explore the intricacies of the EGCG-p53 binding mechanism. A triplicate microsecond MD simulation for each system is initially performed to capture diverse p53 NTD conformations. From the start, the most populated cluster of the second run (R2-1) stands out due to a unique opening between Trp23 and Trp53. During MD simulations, this conformation allows EGCG to sustain a high level of stability and affinity while interacting with both regions of interest and deepening the binding pocket. Structural analysis emphasizes the significance of pyrogallol motifs in EGCG binding. Therefore, the conformational shift in this gap is pivotal, enabling EGCG to impede p53 interactions and manifest its anticancer properties. These findings enhance the present comprehension of the anticancer properties of green tea polyphenols and pave the way for future therapeutic developments.

Ultrasound enhanced diffusion in hydrogels: An experimental and non-equilibrium molecular dynamics study

Focused ultrasound has experimentally been found to enhance the diffusion of nanoparticles; our aim with this work is to study this effect closer using both experiments and non-equilibrium molecular dynamics. Measurements from single particle tracking of 40 nm polystyrene nanoparticles in an agarose hydrogel with and without focused ultrasound are presented and compared with a previous experimental study using 100 nm polystyrene nanoparticles. In both cases, we observed an increase in the mean square displacement during focused ultrasound treatment. We developed a coarse-grained non-equilibrium molecular dynamics model with an implicit solvent to investigate the increase in the mean square displacement and its frequency and amplitude dependencies. This model consists of polymer fibers and two sizes of nanoparticles, and the effect of the focused ultrasound was modeled as an external oscillating force field. A comparison between the simulation and experimental results shows similar mean square displacement trends, suggesting that the particle velocity is a significant contributor to the observed ultrasound-enhanced mean square displacement. The resulting diffusion coefficients from the model are compared to the diffusion equation for a two-time continuous time random walk. The model is found to have the same frequency dependency. At lower particle velocity amplitude values, the model has a quadratic relation with the particle velocity amplitude as described by the two-time continuous time random walk derived diffusion equation, but at higher amplitudes, the model deviates, and its diffusion coefficient reaches the non-hindered diffusion coefficient. This observation suggests that at higher ultrasound intensities in hydrogels, the non-hindered diffusion coefficient can be used.

Reprogramming the Transition States to Enhance C-N Cleavage Efficiency of Rhodococcus opacusl-Amino Acid Oxidase

l-Amino acid oxidase (LAAO) is an important biocatalyst used for synthesizing α-keto acids. LAAO from Rhodococcus opacus (RoLAAO) has a broad substrate spectrum; however, its low total turnover number limits its industrial use. To overcome this, we aimed to employ crystal structure-guided density functional theory calculations and molecular dynamic simulations to investigate the catalytic mechanism. Two key steps were identified: S → [TS1] in step 1 and Int1 → [TS2] in step 2. We reprogrammed the transition states [TS1] and [TS2] to reduce the identified energy barrier and obtain a RoLAAO variant capable of catalyzing 19 kinds of l-amino acids to the corresponding high-value α-keto acids with a high total turnover number, yield (≥95.1 g/L), conversion rate (≥95%), and space-time yields ≥142.7 g/L/d in 12-24 h, in a 5 L reactor. Our results indicated the promising potential of the developed RoLAAO variant for use in the industrial production of α-keto acids while providing a potential catalytic-mechanism-guided protein design strategy to achieve the desired physical and catalytic properties of enzymes.

Systems engineering of Escherichia coli for high-level glutarate production from glucose

Glutarate is a key monomer in polyester and polyamide production. The low efficiency of the current biosynthetic pathways hampers its production by microbial cell factories. Herein, through metabolic simulation, a lysine-overproducing E. coli strain Lys5 is engineered, achieving titer, yield, and productivity of 195.9 g/L, 0.67 g/g glucose, and 5.4 g/L·h, respectively. Subsequently, the pathway involving aromatic aldehyde synthase, monoamine oxidase, and aldehyde dehydrogenase (AMA pathway) is introduced into E. coli Lys5 to produce glutarate from glucose. To enhance the pathway's efficiency, rational mutagenesis on the aldehyde dehydrogenase is performed, resulting in the development of variant Mu5 with a 50-fold increase in catalytic efficiency. Finally, a glutarate tolerance gene cbpA is identified and genomically overexpressed to enhance glutarate productivity. With enzyme expression optimization, the glutarate titer, yield, and productivity of E. coli AMA06 reach 88.4 g/L, 0.42 g/g glucose, and 1.8 g/L·h, respectively. These findings hold implications for improving glutarate biosynthesis efficiency in microbial cell factories.

Molecular Dynamics as a Tool for Virtual Ligand Screening

Rational drug design is essential for new drugs to emerge, especially when the structure of a target protein or nucleic acid is known. To that purpose, high-throughput virtual ligand screening campaigns aim at discovering computationally new binding molecules or fragments to modulate particular biomolecular interactions or biological activities, related to a disease process. The structure-based virtual ligand screening process primarily relies on docking methods which allow predicting the binding of a molecule to a biological target structure with a correct conformation and the best possible affinity. The docking method itself is not sufficient as it suffers from several and crucial limitations (lack of full protein flexibility information, no solvation and ion effects, poor scoring functions, and unreliable molecular affinity estimation).At the interface of computer techniques and drug discovery, molecular dynamics (MD) allows introducing protein flexibility before or after a docking protocol, refining the structure of protein-drug complexes in the presence of water, ions, and even in membrane-like environments, describing more precisely the temporal evolution of the biological complex and ranking these complexes with more accurate binding energy calculations. In this chapter, we describe the up-to-date MD, which plays the role of supporting tools in the virtual ligand screening (VS) process.Without a doubt, using docking in combination with MD is an attractive approach in structure-based drug discovery protocols nowadays. It has proved its efficiency through many examples in the literature and is a powerful method to significantly reduce the amount of required wet experimentations (Tarcsay et al, J Chem Inf Model 53:2990-2999, 2013; Barakat et al, PLoS One 7:e51329, 2012; De Vivo et al, J Med Chem 59:4035-4061, 2016; Durrant, McCammon, BMC Biol 9:71-79, 2011; Galeazzi, Curr Comput Aided Drug Des 5:225-240, 2009; Hospital et al, Adv Appl Bioinforma Chem 8:37-47, 2015; Jiang et al, Molecules 20:12769-12786, 2015; Kundu et al, J Mol Graph Model 61:160-174, 2015; Mirza et al, J Mol Graph Model 66:99-107, 2016; Moroy et al, Future Med Chem 7:2317-2331, 2015; Naresh et al, J Mol Graph Model 61:272-280, 2015; Nichols et al, J Chem Inf Model 51:1439-1446, 2011; Nichols et al, Methods Mol Biol 819:93-103, 2012; Okimoto et al, PLoS Comput Biol 5:e1000528, 2009; Rodriguez-Bussey et al, Biopolymers 105:35-42, 2016; Sliwoski et al, Pharmacol Rev 66:334-395, 2014).

Coenzyme-binding pathway on glutamate dehydrogenase suggested from multiple-binding sites visualized by cryo-electron microscopy

The structure of hexameric glutamate dehydrogenase (GDH) in the presence of the coenzyme nicotinamide adenine dinucleotide phosphate (NADP) was visualized using cryogenic transmission electron microscopy to investigate the ligand-binding pathways to the active site of the enzyme. Each subunit of GDH comprises one hexamer-forming core domain and one nucleotide-binding domain (NAD domain), which spontaneously opens and closes the active-site cleft situated between the two domains. In the presence of NADP, the potential map of GDH hexamer, assuming D3 symmetry, was determined at a resolution of 2.4 Å, but the NAD domain was blurred due to the conformational variety. After focused classification with respect to the NAD domain, the potential maps interpreted as NADP molecules appeared at five different sites in the active-site cleft. The subunits associated with NADP molecules were close to one of the four metastable conformations in the unliganded state. Three of the five binding sites suggested a pathway of NADP molecules to approach the active-site cleft for initiating the enzymatic reaction. The other two binding modes may rarely appear in the presence of glutamate, as demonstrated by the reaction kinetics. Based on the visualized structures and the results from the enzymatic kinetics, we discussed the binding modes of NADP to GDH in the absence and presence of glutamate.

Increased intracellular diffusivity of macromolecules within a mammalian cell by low-intensity pulsed ultrasound

Whilst a number of studies have demonstrated that low-intensity pulsed ultrasound (LIPUS) is a promising therapeutic ultrasound technique that can be used for delivering mild mechanical stimuli to target tissue non-invasively, the underlying biophysical mechanisms still remain unclear. Most mechanism studies have focused explicitly on the effects of LIPUS on the cell membrane and mechanosensitive receptors. In the present study, we propose an additional mechanism by which LIPUS propagation through living cells may directly impact intracellular dynamics, particularly the diffusion transport of biomolecules. To support our hypothesis, human epithelial-like cells (SaOS-2 and HeLa) seeded on a confocal dish placed on a microscope stage were exposed to LIPUS with various exposure conditions (ultrasound frequencies of 0.5, 1 and 3 MHz, peak acoustic pressure of 200 and 400 kPa, a pulse repetition frequency of 1 kHz and a 20 % duty cycle), and the diffusivities of various sizes of biomolecules in the cytoplasm area were measured using fluorescence recovery after photobleaching (FRAP). Furthermore, giant unilamellar vesicles (GUVs) filled with macromolecules were used to examine the physical causal relationship between LIPUS and molecular diffusion changes. Nucleocytoplasmic transport coefficients were also measured by modified FRAP that bleaches the whole cell nuclear region. Extracellular signal-regulated kinases (ERK) activity (the phosphorylation dynamics) was monitored using fluorescence resonance energy transfer (FRET) microscopy. All the measurements were taken during, before and after the LIPUS exposure. Our experimental results clearly showed that the diffusion coefficients of macromolecules within the cell increased with acoustic pressure by 12.1 to 33.5 % during the sonication, and the increments were proportional to their molecular sizes regardless of the ultrasound frequency used. This observation in living cells was consistent with the GUVs exposed to the LIPUS, which indicated that the diffusivity increase was a passive physical response to the acoustic energy of LIPUS. Under the 1 MHz LIPUS exposure with 400 kPa, the passive nucleocytoplasmic transport of enhanced green fluorescent protein (EGFP) was accelerated by 21.4 %. With the same LIPUS exposure condition, both the diffusivity and phosphorylation of ERK induced by EGF treatment were significantly elevated simultaneously, which implied that LIPUS could also modify the kinase kinetics in the signal transduction process. Taken together, this study is the first attempt to uncover the physical link between LIPUS and the dynamics of intracellular macromolecules and related biological processes that LIPUS can possibly increase the diffusivity of intracellular macromolecules, leading to the changes in the basic cellular processes: passive nucleocytoplasmic transport and ERK. Our findings can provide a novel perspective that the mechanotransduction process that the intracellular region, in addition to the cell membrane, can convert the acoustic stimuli of LIPUS to biochemical signals.

A Hybrid Structure-Based Machine Learning Approach for Predicting Kinase Inhibition by Small Molecules

Kinases have been the focus of drug discovery programs for three decades leading to over 70 therapeutic kinase inhibitors and biophysical affinity measurements for over 130,000 kinase-compound pairs. Nonetheless, the precise target spectrum for many kinases remains only partly understood. In this study, we describe a computational approach to unlocking qualitative and quantitative kinome-wide binding measurements for structure-based machine learning. Our study has three components: (i) a Kinase Inhibitor Complex (KinCo) data set comprising in silico predicted kinase structures paired with experimental binding constants, (ii) a machine learning loss function that integrates qualitative and quantitative data for model training, and (iii) a structure-based machine learning model trained on KinCo. We show that our approach outperforms methods trained on crystal structures alone in predicting binary and quantitative kinase-compound interaction affinities; relative to structure-free methods, our approach also captures known kinase biochemistry and more successfully generalizes to distant kinase sequences and compound scaffolds.

Crosslinking Rapidly Cured Epoxy Resin Thermosets: Experimental and Computational Modeling and Simulation Study

The power of computational modeling and simulation for establishing clear links between materials' intrinsic properties and their atomic structure has more and more increased the demand for reliable and reproducible protocols. Despite this increased demand, no one approach can provide reliable and reproducible outcomes to predict the properties of novel materials, particularly rapidly cured epoxy-resins with additives. This study introduces the first computational modeling and simulation protocol for crosslinking rapidly cured epoxy resin thermosets based on solvate ionic liquid (SIL). The protocol combines several modeling approaches, including quantum mechanics (QMs) and molecular dynamics (MDs). Furthermore, it insightfully provides a wide range of thermo-mechanical, chemical, and mechano-chemical properties, which agree with experimental data.

Advancing the field of computational drug design using multicanonical molecular dynamics-based dynamic docking

Multicanonical molecular dynamics (McMD)-based dynamic docking is a powerful tool to not only predict the native binding configuration between two flexible molecules, but it can also be used to accurately simulate the binding/unbinding pathway. Furthermore, it can also predict alternative binding sites, including allosteric ones, by employing an exhaustive sampling approach. Since McMD-based dynamic docking accurately samples binding/unbinding events, it can thus be used to determine the molecular mechanism of binding between two molecules. We developed the McMD-based dynamic docking methodology based on the powerful, but woefully underutilized McMD algorithm, combined with a toolset to perform the docking and to analyze the results. Here, we showcase three of our recent works, where we have applied McMD-based dynamic docking to advance the field of computational drug design. In the first case, we applied our method to perform an exhaustive search between Hsp90 and one of its inhibitors to successfully predict the native binding configuration in its binding site, as we refined our analysis methods. For our second case, we performed an exhaustive search of two medium-sized ligands and Bcl-xL, which has a cryptic binding site that differs greatly between the apo and holo structures. Finally, we performed a dynamic docking simulation between a membrane-embedded GPCR molecule and a high affinity ligand that binds deep within its receptor's pocket. These advanced simulations showcase the power that the McMD-based dynamic docking method has, and provide a glimpse of the potential our methodology has to unravel and solve the medical and biophysical issues in the modern world.

Supplementary Information

The online version contains supplementary material available at 10.1007/s12551-022-01010-z.

Ligand-induced shifts in conformational ensembles that describe transcriptional activation

Nuclear receptors function as ligand-regulated transcription factors whose ability to regulate diverse physiological processes is closely linked with conformational changes induced upon ligand binding. Understanding how conformational populations of nuclear receptors are shifted by various ligands could illuminate strategies for the design of synthetic modulators to regulate specific transcriptional programs. Here, we investigate ligand-induced conformational changes using a reconstructed, ancestral nuclear receptor. By making substitutions at a key position, we engineer receptor variants with altered ligand specificities. We combine cellular and biophysical experiments to characterize transcriptional activity, as well as elucidate mechanisms underlying altered transcription in receptor variants. We then use atomistic molecular dynamics (MD) simulations with enhanced sampling to generate ensembles of wildtype and engineered receptors in combination with multiple ligands, followed by conformational analysis and correlation of MD-based predictions with functional ligand profiles. We determine that conformational ensembles accurately describe ligand responses based on observed population shifts. These studies provide a platform which will allow structural characterization of physiologically-relevant conformational ensembles, as well as provide the ability to design and predict transcriptional responses in novel ligands.

Computationally Elucidating the Binding Kinetics for Different AChE Inhibitors to Access the Rationale for Improving the Drug Efficacy

Traditional drug discovery is based on a binding affinity (thermodynamics)-driven paradigm. Numerous examples, however, demonstrated that drug efficacy does not always depend only on binding affinity but positively correlates with binding kinetics, that is, the dissociation rate constant (k_off). Binding free energy landscape (BFEL) constructor is a computational binding kinetics prediction method, previously developed by us, that estimates the binding kinetics for ligand-protein based on their constructed binding free energy landscape, but it also reveals the detailed molecular mechanism of the binding event, hence, providing the position of transition states at the molecular level to modify/improve the binding kinetics. Acetylcholinesterase (AChE) is a well-known Alzheimer's disease (AD) target for which there is still not an ideal drug on the market. Therefore, to improve the drug design strategy for AD, the binding kinetics and binding molecular mechanisms of the four inhibitors of AChE, that is, E2020 (Aricept), HupA, Rivastigmine, and Galantamine, were studied. Also, the differentiation of the binding kinetics between mAChE and TcAChE was studied to evaluate the sensitiveness of BFEL constructor. The flexibility of molecules has a noticeable effect on the nature of BFEL. To the same target, flexible molecules (i.e., E2020 and Rivastigmine) which contain more rotatable bonds tend to have more complicated BFELs reflecting more complicated molecular action mechanisms than the rigid ones (i.e., HupA and Galantamine), which therefore could be more challenging to be optimized. The binding kinetics is highly dependent on the structure of the molecules, such as the length and the functional groups. Therefore, E2020 presents better binding kinetic and thermodynamic properties with either TcAChE or mAChE. Therefore, it is the most promising lead drug for binding kinetics-based drug design. In addition, the binding kinetics of a drug may present different values in the proteins of different organisms because the residue compositions of the binding gorges of the targets are variant, that is, E2020 shows lower binding affinity and association energy barrier in binding with mAChE than TcAChE. However, HupA presents a better binding property with TcAChE than mAChE.

Substrate Tumbling in a Chemisorbed Diastereomeric α-Ketoester/1-(1-Naphthyl)ethylamine Complex

Scanning tunneling microscopy (STM) data for α-ketoester/1-(1-naphthyl)ethylamine complexes on Pt(111) reveal a tumbling motion that couples two neighboring binding states. The interconversion, resulting in prochiral inversion of the α-ketoester, occurs in single complexes without breaking them apart. This is a surprising observation because the overall motion requires rotation of the α-ketoester away from the surface without branching exclusively into diffusion away from the complex or desorption. The multi-step interconversion is rationalized in terms of sequences of bound states that combine transient H-bond interactions with the chiral molecule and weakened adsorption interactions with the metal. The observation of tumbling in single long-lived complexes is of relevance to self-assembly and directed molecular motion on surfaces, to ligand-controlled surface reactions, and most directly to stereocontrol in asymmetric heterogeneous catalysis.

On the force field optimisation of β -lactam cores using the force field Toolkit

When employing molecular dynamics (MD) simulations for computer-aided drug design, the quality of the used force fields is highly important. Here we present reparametrisations of the force fields for the core molecules from 9 differentβ -lactam classes, for which we utilized the force field Toolkit and Gaussian calculations. We focus on the parametrisation of the dihedral angles, with the goal of reproducing the optimised quantum geometry in MD simulations. Parameters taken from CGenFF turn out to be a good initial guess for the multiplicity of each dihedral angle, but the key to a successful parametrisation is found to lie in the phase shifts. Based on the optimised quantum geometry, we come up with a strategy for predicting the phase shifts prior to the dihedral potential fitting. This allows us to successfully parameterise 8 out of the 11 molecules studied here, while the remaining 3 molecules can also be parameterised with small adjustments. Our work highlights the importance of predicting the dihedral phase shifts in the ligand parametrisation protocol, and provides a simple yet valuable strategy for improving the process of parameterising force fields of drug-like molecules.

In Silico Modeling and Scoring of PROTAC-Mediated Ternary Complex Poses

Proteolysis targeting chimeras (PROTACs) are molecules that induce protein degradation via formation of ternary complexes between an E3 ubiquitin ligase and a target protein. The rational design of PROTACs requires accurate knowledge of the native configuration of the PROTAC-induced ternary complex. This study demonstrates that native and non-native ternary complex poses can be distinguished based on the pose occupancy time in MD, where native poses exhibit longer occupancy times at both room and higher temperatures. Candidate poses are generated by MD sampling and pre-ranked by classic MM/GBSA. A specific heating scheme is then applied to accelerate ternary pose departure, with the pose occupancy time and fraction being measured. This scoring identifies the native pose in all systems tested. Its success is partially attributed to the dynamic nature of pose departure analyses, which accounts for entropic effects typically neglected in the faster static scoring methods, while entropy plays a greater role in protein-protein than in protein-ligand systems.

From Data to Knowledge: Systematic Review of Tools for Automatic Analysis of Molecular Dynamics Output

An increasing number of crystal structures available on one side, and the boost of computational power available for computer-aided drug design tasks on the other, have caused that the structure-based drug design tools are intensively used in the drug development pipelines. Docking and molecular dynamics simulations, key representatives of the structure-based approaches, provide detailed information about the potential interaction of a ligand with a target receptor. However, at the same time, they require a three-dimensional structure of a protein and a relatively high amount of computational resources. Nowadays, as both docking and molecular dynamics are much more extensively used, the amount of data output from these procedures is also growing. Therefore, there are also more and more approaches that facilitate the analysis and interpretation of the results of structure-based tools. In this review, we will comprehensively summarize approaches for handling molecular dynamics simulations output. It will cover both statistical and machine-learning-based tools, as well as various forms of depiction of molecular dynamics output.

Recent Developments in Positron Emission Tomography Tracers for Proteinopathies Imaging in Dementia

An early detection and intervention for dementia represent tremendous unmet clinical needs and priorities in society. A shared feature of neurodegenerative diseases causing dementia is the abnormal accumulation and spreading of pathological protein aggregates, which affect the selective vulnerable circuit in a disease-specific pattern. The advancement in positron emission tomography (PET) biomarkers has accelerated the understanding of the disease mechanism and development of therapeutics for Alzheimer's disease and Parkinson's disease. The clinical utility of amyloid-β PET and the clinical validity of tau PET as diagnostic biomarker for Alzheimer's disease continuum have been demonstrated. The inclusion of biomarkers in the diagnostic criteria has introduced a paradigm shift that facilitated the early and differential disease diagnosis and impacted on the clinical management. Application of disease-modifying therapy likely requires screening of patients with molecular evidence of pathological accumulation and monitoring of treatment effect assisted with biomarkers. There is currently still a gap in specific 4-repeat tau imaging probes for 4-repeat tauopathies and α-synuclein imaging probes for Parkinson's disease and dementia with Lewy body. In this review, we focused on recent development in molecular imaging biomarkers for assisting the early diagnosis of proteinopathies (i.e., amyloid-β, tau, and α-synuclein) in dementia and discussed future perspectives.

N-Terminal-Driven Binding Mechanism of an Antigen Peptide to Human Leukocyte Antigen-A*2402 Elucidated by Multicanonical Molecular Dynamic-Based Dynamic Docking and Path Sampling Simulations

We have applied our advanced multicanonical molecular dynamics (McMD)-based dynamic docking methodology to investigate the binding mechanism of an HIV-1 Nef protein epitope to the Asian-dominant allele human leukocyte antigen (HLA)-A*2402. Even though pMHC complex formation [between a Major histocompatibility complex (MHC) class I molecule, which is encoded by an HLA allele, and an antigen peptide] is one of the fundamental processes of the adaptive human immune response, its binding mechanism has not yet been well studied, partially due to the high allelic variation of HLAs in the population. We have used our developed McMD-based dynamic docking method and have successfully reproduced the native complex structure, which is located near the free energy global minimum. Subsequent path sampling MD simulations elucidated the atomic details of the binding process and indicated that the peptide binding is initially driven by the highly positively charged N-terminus of the peptide that is attracted to the various negatively charged residues on the MHC molecule's surface. Upon nearing the pocket, the second tyrosine residue of the peptide anchors the peptide by strongly binding to the B-site of the MHC molecule via hydrophobic driven interactions, resulting in a very strong bound complex structure. Our methodology can be effectively used to predict the bound complex structures between MHC molecules and their antigens to study their binding mechanism in close detail, which would help with the development of new vaccines against cancers, as well as viral infections such as HIV and COVID-19.

Accurate Binding Configuration Prediction of a G-Protein-Coupled Receptor to Its Antagonist Using Multicanonical Molecular Dynamics-Based Dynamic Docking

We have performed dynamic docking between a prototypic G-protein-coupled receptor (GPCR) system, the β₂-adrenergic receptor, and its antagonist, alprenolol, using one of the enhanced conformation sampling methods, multicanonical molecular dynamics (McMD), which does not rely on any prior knowledge for the definition of the reaction coordinate. Although we have previously applied our McMD-based dynamic docking protocol to various globular protein systems, its application to GPCR systems would be difficult because of their complicated design, which include a lipid bilayer, and because of the difficulty in sampling the configurational space of a binding site that exists deep inside the GPCR. Our simulations sampled a wide array of ligand-bound and ligand-unbound structures, and we measured 427 binding events during our 48 μs production run. Analysis of the ensemble revealed several stable and meta-stable structures, where the most stable structure at the global free energy minimum matches the experimental one. Additional canonical MD simulations were used for refinement and validation of the structures, revealing that most of the intermediates are sufficiently stable to trap the ligand in these intermediary states and furthermore validated our prediction results. Given the difficulty in reaching the orthosteric binding site, chemical optimization of the compound for the second ranking configuration, which binds near the pocket's entrance, might lead to a high-affinity allosteric inhibitor. Accordingly, we show that the application of our methodology can be used to provide crucial insights for the rational design of drugs that target GPCRs.

Probing the Suitability of Different Ca2+ Parameters for Long Simulations of Diisopropyl Fluorophosphatase

Organophosphate hydrolases are promising as potential biotherapeutic agents to treat poisoning with pesticides or nerve gases. However, these enzymes often need to be further engineered in order to become useful in practice. One example of such enhancement is the alteration of enantioselectivity of diisopropyl fluorophosphatase (DFPase). Molecular modeling techniques offer a unique opportunity to address this task rationally by providing a physical description of the substrate-binding process. However, DFPase is a metalloenzyme, and correct modeling of metal cations is a challenging task generally coming with a tradeoff between simulation speed and accuracy. Here, we probe several molecular mechanical parameter combinations for their ability to empower long simulations needed to achieve a quantitative description of substrate binding. We demonstrate that a combination of the Amber19sb force field with the recently developed 12-6 Ca2+ models allows us to both correctly model DFPase and obtain new insights into the DFP binding process.

Pepsin-like aspartic proteases (PAPs) as model systems for combining biomolecular simulation with biophysical experiments

Pepsin-like aspartic proteases (PAPs) are a class of aspartic proteases which shares tremendous structural similarity with human pepsin. One of the key structural features of PAPs is the presence of a β-hairpin motif otherwise known as flap. The biological function of the PAPs is highly dependent on the conformational dynamics of the flap region. In apo PAPs, the conformational dynamics of the flap is dominated by the rotational degrees of freedom associated with χ1 and χ2 angles of conserved Tyr (or Phe in some cases). However it is plausible that dihedral order parameters associated with several other residues might play crucial roles in the conformational dynamics of apo PAPs. Due to their size, complexities associated with conformational dynamics and clinical significance (drug targets for malaria, Alzheimer's disease etc.), PAPs provide a challenging testing ground for computational and experimental methods focusing on understanding conformational dynamics and molecular recognition in biomolecules. The opening of the flap region is necessary to accommodate substrate/ligand in the active site of the PAPs. The BIG challenge is to gain atomistic details into how reversible ligand binding/unbinding (molecular recognition) affects the conformational dynamics. Recent reports of kinetics (K i, K d) and thermodynamic parameters (ΔH, TΔS, and ΔG) associated with macro-cyclic ligands bound to BACE1 (belongs to PAP family) provide a perfect challenge (how to deal with big ligands with multiple torsional angles and select optimum order parameters to study reversible ligand binding/unbinding) for computational methods to predict binding free energies and kinetics beyond typical test systems e.g. benzamide-trypsin. In this work, i reviewed several order parameters which were proposed to capture the conformational dynamics and molecular recognition in PAPs. I further highlighted how machine learning methods can be used as order parameters in the context of PAPs. I then proposed some open ideas and challenges in the context of molecular simulation and put forward my case on how biophysical experiments e.g. NMR, time-resolved FRET etc. can be used in conjunction with biomolecular simulation to gain complete atomistic insights into the conformational dynamics of PAPs.