MEPSA: minimum energy pathway analysis for energy landscapes

Metadynamics Study of Lipid-Mediated Antibacterial Toxin Binding to the EmrE Multiefflux Protein

EmrE is a bacterial efflux protein in the small multidrug-resistant (SMR) family present in Escherichia coli. Due to its small size, 110 residues in each dimer subunit, it is an ideal model system to study ligand-protein-membrane interactions. Here in our work, we have calculated the free energy landscape of benzyltrimetylammonium (BTMA) and tetraphenyl phosphonium (TPP) binding to EmrE using the enhanced sampling method-multiple walker metadynamics. We estimate that the free energy of BTMA binding to EmrE is -21.2 ± 3.3 kJ/mol and for TPP is -43.6 ± 3.8 kJ/mol. BTMA passes through two metastable states to reach the binding pocket, while TPP has a more complex binding landscape with four metastable states and one main binding site. Our simulations show that the ligands interact with the membrane lipids at a distance 1 nm away from the binding site which forms a broad local minimum, consistent for both BTMA and TPP. This site can be an alternate entry point for ligands to partition from the membrane into the protein, especially for bulky and/or branched ligands. We also observed the membrane lipid and C-terminal 110HisA form salt-bridge interactions with the helix-1 residue 22LysB. Our free energy estimates and clusters are in close agreement with experimental data and give us an atomistic view of the ligand-protein-lipid interactions. Understanding the binding pathway of these ligands can guide us in future design of ligands that can alter or halt the function of EmrE.

Structural basis of the obligatory exchange mode of human neutral amino acid transporter ASCT2

ASCT2 is an obligate exchanger of neutral amino acids, contributing to cellular amino acid homeostasis. ASCT2 belongs to the same family (SLC1) as Excitatory Amino Acid Transporters (EAATs) that concentrate glutamate in the cytosol. The mechanism that makes ASCT2 an exchanger rather than a concentrator remains enigmatic. Here, we employ cryo-electron microscopy and molecular dynamics simulations to elucidate the structural basis of the exchange mechanism of ASCT2. We establish that ASCT2 binds three Na+ ions per transported substrate and visits a state that likely acts as checkpoint in preventing Na+ ion leakage, both features shared with EAATs. However, in contrast to EAATs, ASCT2 retains one Na+ ion even under Na+-depleted conditions. We demonstrate that ASCT2 cannot undergo the structural transition in TM7 that is essential for the concentrative transport cycle of EAATs. This structural rigidity and the high-affinity Na+ binding site effectively confine ASCT2 to an exchange mode.

Mechanistic insights into the conversion of flavin adenine dinucleotide (FAD) to 8-formyl FAD in formate oxidase: a combined experimental and in-silico study

Formate oxidase (FOx), which contains 8-formyl flavin adenine dinucleotide (FAD), exhibits a distinct advantage in utilizing ambient oxygen molecules for the oxidation of formic acid compared to other glucose-methanol-choline (GMC) oxidoreductase enzymes that contain only the standard FAD cofactor. The FOx-mediated conversion of FAD to 8-formyl FAD results in an approximate 10-fold increase in formate oxidase activity. However, the mechanistic details underlying the autocatalytic formation of 8-formyl FAD are still not well understood, which impedes further utilization of FOx. In this study, we employ molecular dynamics simulation, QM/MM umbrella sampling simulation, enzyme activity assay, site-directed mutagenesis, and spectroscopic analysis to elucidate the oxidation mechanism of FAD to 8-formyl FAD. Our results reveal that a catalytic water molecule, rather than any catalytic amino acids, serves as a general base to deprotonate the C8 methyl group on FAD, thus facilitating the formation of a quinone-methide tautomer intermediate. An oxygen molecule subsequently oxidizes this intermediate, resulting in a C8 methyl hydroperoxide anion that is protonated and dissociated to form OHC-RP and OH-. During the oxidation of FAD to 8-formyl FAD, the energy barrier for the rate-limiting step is calculated to be 22.8 kcal/mol, which corresponds to the required 14-hour transformation time observed experimentally. Further, the elucidated oxidation mechanism reveals that the autocatalytic formation of 8-formyl FAD depends on the proximal arginine and serine residues, R87 and S94, respectively. Enzymatic activity assay validates that the mutation of R87 to lysine reduces the kcat value to 75% of the wild-type, while the mutation to histidine results in a complete loss of activity. Similarly, the mutant S94I also leads to the deactivation of enzyme. This dependency arises because the nucleophilic OH- group and the quinone-methide tautomer intermediate are stabilized through the noncovalent interaction provided by R87 and S94. These findings not only explain the mechanistic details of each reaction step but also clarify the functional role of R87 and S94 during the oxidative maturation of 8-formyl FAD, thereby providing crucial theoretical support for the development of novel flavoenzymes with enhanced redox properties.

Mechanistic variances in NO release: ortho vs. meta isomers of nitrophenol and nitroaniline

The NO release following 266 nm photolysis of ortho and meta isomers of nitrophenol and nitroaniline shows a bimodal translational energy distribution, wherein the slow and fast components originate from dynamics in the S0 and T1 states, respectively. The translational energy distribution profiles for any NO product state show a higher slow-to-fast (s/f) branching ratio for the ortho isomer in comparison with the meta isomer. The observed variation in the s/f branching ratio vis-à-vis the ortho and meta isomers is attributed to the presence of intramolecular hydrogen bonding between the ortho substituent and NO2 moiety, which favours the roaming mechanism.

Computational Analysis of Activation of Dimerized Epidermal Growth Factor Receptor Kinase Using the String Method and Markov State Model

Epidermal growth factor receptor (EGFR) activation is accompanied by dimerization. During the activation of the intracellular kinase domain, two EGFR kinases form an asymmetric dimer, and one side of the dimer (receiver) is activated. Using the string method and Markov state model (MSM), we performed a computational analysis of the structural changes in the activation of the EGFR dimer in this study. The string method reveals the minimum free-energy pathway (MFEP) from the inactive to active structure. The MSM was constructed from numerous trajectories of molecular dynamics simulations around the MFEP, which revealed the free-energy map of structural changes. In the activation of the receiver kinase, the unfolding of the activation loop (A-loop) is followed by the rearrangement of the C-helix, as observed in other kinases. However, unlike other kinases, the free-energy map of EGFR at the asymmetric dimer showed that the active state yielded the highest stability and revealed how interactions at the dimer interface induced receiver activation. As the H-helix of the activator approaches the C-helix of the receiver during activation, the A-loop unfolds. Subsequently, L782 of the receiver enters the pocket between the G- and H-helices of the activator, leading to a rearrangement of the hydrophobic residues around L782 of the receiver, which constitutes a structural rearrangement of the C-helix of the receiver from an outward to an inner position. The MSM analysis revealed long-time scale trajectories via kinetic Monte Carlo.

Can local heating and molecular crowders disintegrate amyloid aggregates?

The present study employs a blend of molecular dynamics simulations and a theoretical model to explore the potential disintegration mechanism of a matured Aβ octamer, aiming to offer a strategy to combat Alzheimer's disease. We investigate local heating and crowding effects on Aβ disintegration by selectively heating key Aβ segments and varying the concentration of sodium dodecyl sulphate (SDS), respectively. Despite initiation of disruption, Aβ aggregates resist complete disintegration during local heating due to rapid thermal energy distribution to the surrounding water. Conversely, although SDS molecules effectively inhibit Aβ aggregation at higher concentration through micelle formation, they fail to completely disintegrate the aggregate due to the exceedingly high energy barrier. To address the sampling challenge posed by the formidable energy barrier, we have performed well-tempered metadynamics simulations. Simulations reveal a multi-step disintegration mechanism for the Aβ octamer, suggesting a probable sequence: octamer → pentamer/hexamer ⇌ tetramer → monomer, with a rate-determining step constituting 45 kJ mol-1 barrier during the octamer to pentamer/hexamer transition. Additionally, we have proposed a novel two-state mean-field model based on Ising spins that offers an insight into the kinetics of the Aβ growth process and external perturbation effects on disintegration. Thus, the current simulation study, coupled with the newly introduced mean-field model, offers an insight into the detailed mechanisms underlying the Aβ aggregation process, guiding potential strategies for effective disintegration of Aβ aggregates.

On-Surface Synthesis of Non-Benzenoid Nanographenes Embedding Azulene and Stone-Wales Topologies

The incorporation of non-benzenoid motifs in graphene nanostructures significantly impacts their properties, making them attractive for applications in carbon-based electronics. However, understanding how specific non-benzenoid structures influence their properties remains limited, and further investigations are needed to fully comprehend their implications. Here, we report an on-surface synthetic strategy toward fabricating non-benzenoid nanographenes containing different combinations of pentagonal and heptagonal rings. Their structure and electronic properties were investigated via scanning tunneling microscopy and spectroscopy, complemented by computational investigations. After thermal activation of the precursor P on the Au(111) surface, we detected two major nanographene products. Nanographene Aa-a embeds two azulene units formed through oxidative ring-closure of methyl substituents, while Aa-s contains one azulene unit and one Stone-Wales defect, formed by the combination of oxidative ring-closure and skeletal ring-rearrangement reactions. Aa-a exhibits an antiferromagnetic ground state with the highest magnetic exchange coupling reported up to date for a non-benzenoid containing nanographene, coexisting with side-products with closed shell configurations resulted from the combination of ring-closure and ring-rearragement reactions (Ba-a , Ba-s , Bs-a and Bs-s ). Our results provide insights into the single gold atom assisted synthesis of novel NGs containing non-benzenoid motifs and their tailored electronic/magnetic properties.

An alpha-helical lid guides the target DNA toward catalysis in CRISPR-Cas12a

CRISPR-Cas12a is a powerful RNA-guided genome-editing system that generates double-strand DNA breaks using its single RuvC nuclease domain by a sequential mechanism in which initial cleavage of the non-target strand is followed by target strand cleavage. How the spatially distant DNA target strand traverses toward the RuvC catalytic core is presently not understood. Here, continuous tens of microsecond-long molecular dynamics and free-energy simulations reveal that an α-helical lid, located within the RuvC domain, plays a pivotal role in the traversal of the DNA target strand by anchoring the crRNA:target strand duplex and guiding the target strand toward the RuvC core, as also corroborated by DNA cleavage experiments. In this mechanism, the REC2 domain pushes the crRNA:target strand duplex toward the core of the enzyme, while the Nuc domain aids the bending and accommodation of the target strand within the RuvC core by bending inward. Understanding of this critical process underlying Cas12a activity will enrich fundamental knowledge and facilitate further engineering strategies for genome editing.

Probing the pH Sensitivity of OprM: Insights into Metastable States and Semi-Open Conformation

Efflux pumps are specialized transport proteins that play a key role in the bacterial defense against a wide spectrum of antibiotics. Hence, understanding the biophysical mechanism associated with this complex system of drug expulsion becomes crucial. This work deals with some vital aspects of the outer membrane factor (OMF) of MexAB-OprM. After being passed through MexB and MexA, efflux substrates have to go through OprM for their final judgment. Thus, it is very important to understand the periplasmic pore opening mechanism and the associated biophysical changes during this process. Our study captures a detailed analysis of the pore opening mechanism involving OprM. With powerful molecular dynamics (MD) techniques such as well-tempered metadynamics, the presence of metastable states in between open and closed states was confirmed. Also, upon mutating R376, the energy barrier for the conversion of the close to open conformation decreases, indicating an important role played by the residue. Further, constant pH MD was performed to capture the effect of pH in both conformations. OprM exhibits distinct conformational states at pH values greater than 5.5 and lower than 5.5, suggesting its pH-responsive characteristics. Overall, our study elucidates a crucial undertaking toward discovering potential inhibitors for MexAB-OprM efflux pumps.

Insight into the stacking effect on shifted patterns of bilayer phosphorene: a comprehensive first-principles study

It is crucial to deeply understand how the interlayer interaction acts on controlling the structural and electronic properties of shifted patterns of bilayer phosphorene. A comprehensive first-principles study on the bilayer phosphorene through relative translation along different directions has revealed that there is a direct correlation between the potential energy surface and the interlayer equilibrium distance. The shorter the interlayer distance, the lower the potential energy surface. The shifted patterns with the most stable state, the metastable state, and the transition state (with energy barrier of ∼1.3 meV/atom) were found associated with the AB, the Aδ, and the TS stacking configurations, respectively. The high energy barriers, on the other hand, are ∼9.3 meV/atom at the AA stacking configuration along the zigzag pathway, ∼5.3 meV/atom at the AB' stacking configuration along the armchair pathway, and ∼11.2 meV/atom at the AA' stacking configuration along the diagonal pathway, respectively. The character of electronic bandgap with respect to the shifting shows an anisotropic behavior (with the value of 0.69-1.22 eV). A transition from the indirect to the direct bandgap occurs under the shifting, implying a tunable bandgap by stacking engineering. Furthermore, the orbital hybridization at the interfacial region induces a redistribution of the net charge (∼0.002-0.011e) associated with the relative shifting between layers, leading to a strong polarization with stripe-like electron depletion near the lone pairs and accumulation in the middle of the interfacial region. It is expected that such interesting findings will provide a fundamental reference to deeply understand and analyze the complex local structural and electronic properties of twisted bilayer phosphorene and will make the shifted patterns of bilayer phosphorene promising for nanoelectronics as versatile shiftronics materials.

The role of dynamics in heterogeneous catalysis: Surface diffusivity and N2 decomposition on Fe(111)

Dynamics has long been recognized to play an important role in heterogeneous catalytic processes. However, until recently, it has been impossible to study their dynamical behavior at industry-relevant temperatures. Using a combination of machine learning potentials and advanced simulation techniques, we investigate the cleavage of the N[Formula: see text] triple bond on the Fe(111) surface. We find that at low temperatures our results agree with the well-established picture. However, if we increase the temperature to reach operando conditions, the surface undergoes a global dynamical change and the step structure of the Fe(111) surface is destabilized. The catalytic sites, traditionally associated with this surface, appear and disappear continuously. Our simulations illuminate the danger of extrapolating low-temperature results to operando conditions and indicate that the catalytic activity can only be inferred from calculations that take dynamics fully into account. More than that, they show that it is the transition to this highly fluctuating interfacial environment that drives the catalytic process.

The H163A mutation unravels an oxidized conformation of the SARS-CoV-2 main protease

The main protease of SARS-CoV-2 (Mpro) is an important target for developing COVID-19 therapeutics. Recent work has highlighted Mpro's susceptibility to undergo redox-associated conformational changes in response to cellular and immune-system-induced oxidation. Despite structural evidence indicating large-scale rearrangements upon oxidation, the mechanisms of conformational change and its functional consequences are poorly understood. Here, we present the crystal structure of an Mpro point mutant (H163A) that shows an oxidized conformation with the catalytic cysteine in a disulfide bond. We hypothesize that Mpro adopts this conformation under oxidative stress to protect against over-oxidation. Our metadynamics simulations illustrate a potential mechanism by which H163 modulates this transition and suggest that this equilibrium exists in the wild type enzyme. We show that other point mutations also significantly shift the equilibrium towards this state by altering conformational free energies. Unique avenues of SARS-CoV-2 research can be explored by understanding how H163 modulates this equilibrium.

Photodegradation of Flutamide and Halogen Derivatives of Nitrobenzotrifluoride: The NO Release Channel

The photodegradation of the nonsteroidal antiandrogen drug flutamide has been long linked to the photoisomerization involving the nitro group. In this work, the dynamics of NO photoelimination upon photolysis at 266 nm of flutamide, nitrobenzotrifluoride, and its halogen derivatives were investigated. Similar to nitrobenzene and its derivatives, a bimodal translational energy distribution was observed for the NO photofragment indicating the presence of two distinct elimination channels resulting in slow and fast components. The trends in the slow/fast branching ratio show that halogen substitution at the para position increases the triplet state yield due to the internal heavy-atom effect leading to enhancement of the fast component. Furthermore, the topology of the triplet state potential energy surface showed that the minimum energy path favors the oxaziridine ring-type intermediate over the NO2 roaming mechanism in all five molecules investigated. The steric interaction between the NO2 group and the CF3 group, which are placed in the ortho position, lowers the barrier for the formation of the oxaziridine transition state compared to that of nitrobenzene.

Beyond ManifoldEM: geometric relationships between manifold embeddings of a continuum of 3D molecular structures and their 2D projections

ManifoldEM is an established method of geometric machine learning developed to extract information on conformational motions of molecules from their projections obtained by cryogenic electron microscopy (cryo-EM). In a previous work, in-depth analysis of the properties of manifolds obtained for simulated ground-truth data from molecules exhibiting domain motions has led to improvements of this method, as demonstrated in selected applications of single-particle cryo-EM. In the present work this analysis has been extended to investigate the properties of manifolds constructed by embedding data from synthetic models represented by atomic coordinates in motion, or three-dimensional density maps from biophysical experiments other than single-particle cryo-EM, with extensions to cryo-electron tomography and single-particle imaging with a X-ray free-electron laser. Our theoretical analysis revealed interesting relationships between all these manifolds, which can be exploited in future work.

Modulating the Roaming Dynamics for the NO Release in ortho-Nitrobenzenes

The dynamics of NO release upon photodissociation of nitroaromatic compounds is dependent on the nature of the interaction between the NO₂ group and substituent in the ortho position. A bimodal (slow and fast) translational energy distribution of the NO photofragment indicates the presence of two distinct NO elimination channels. The slow-to-fast branching ratio for the NO release is regulated by the hydrogen bonding ability of the ortho substituent and follows the order [OH > NH₂ > CH₃ > OCH₃], indicating that the intramolecular hydrogen bonding plays a pivotal role in NO release dynamics. Further, the topology of the triplet state potential energy surface acts as a doorway to the dissociation pathway switching between the roaming and nonroaming mechanisms, with hydrogen bonding substituents (OH and NH₂) favoring the roaming mechanism.

Insights into the Degradation of Polymer-Drug Conjugates by an Overexpressed Enzyme in Cancer Cells

Intensive efforts have been made to provide better treatments to cancer patients. Currently, nanoparticle-based drug delivery systems have gained propulsion, as they can overcome the drawbacks of free drugs. However, drug stability inside the nanocapsule must be ensured to prevent burst release. To overcome this, drugs conjugated to amphiphilic copolymers, assembled into nanoparticles, can provide a sustained release if endogenously degraded. Thus, we have designed and assessed the drug release viability of polymer-drug conjugates by the human Carboxylesterase 2, for a targeted drug activation. We performed molecular dynamics simulations applying a quantum mechanics/molecular mechanics potential to study the degradation profiles of 30 designed conjugates, where six were predicted to be hydrolyzed by this enzyme. We further analyzed the enzyme-substrate environment to delve into what structural features may lead to successful hydrolysis. These findings contribute to the development of new medicines ensuring effective cancer treatments with fewer side effects.

Mechanistic Modeling of Lys745 Sulfonylation in EGFR C797S Reveals Chemical Determinants for Inhibitor Activity and Discriminates Reversible from Irreversible Agents

Targeted covalent inhibitors hold promise for drug discovery, particularly for kinases. Targeting the catalytic lysine of epidermal growth factor receptor (EGFR) has attracted attention as a new strategy to overcome resistance due to the emergence of C797S mutation. Sulfonyl fluoride derivatives able to inhibit EGFR^{L858R/T790M/C797S} by sulfonylation of Lys745 have been reported. However, atomistic details of this process are still poorly understood. Here, we describe the mechanism of inhibition of an innovative class of compounds that covalently engage the catalytic lysine of EGFR, through a sulfur(VI) fluoride exchange (SuFEx) process, with the help of hybrid quantum mechanics/molecular mechanics (QM/MM) and path collective variables (PCVs) approaches. Our simulations identify the chemical determinants accounting for the irreversible activity of agents targeting Lys745 and provide hints for the further optimization of sulfonyl fluoride agents.

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

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

Force-Field-Dependent DNA Breathing Dynamics: A Case Study of Hoogsteen Base Pairing in A6-DNA

The Hoogsteen (HG) base pairing conformation, commonly observed in damaged and mutated DNA helices, facilitates DNA repair and DNA recognition. The free energy difference between HG and Watson-Crick (WC) base pairs has been computed in previous studies. However, the mechanism of the conformational transition is not well understood. A detailed understanding of the process of WC to HG base pair transition can provide a deeper understanding of DNA repair and recognition. In an earlier study, we explored the free energy landscape for this process using extensive computer simulation with the CHARMM36 force field. In this work, we study the impact of force field models in describing the WC to HG base pairing transition using meta-eABF enhanced sampling, quasi-harmonic entropy calculation, and nonbonded energy analysis. The secondary structures of both base pairing forms and the topology of the free energy landscapes were consistent over different force field models, although the relative free energy, entropy, and the interaction energies tend to vary. The relative stability of the WC and HG conformations is dictated by a delicate balance between the enthalpic stabilization and the reduced entropy of the structurally rigid HG structure. These findings highlight the impact that subtleties in force field models can have on accurately modeling DNA base pair dynamics and should stimulate further computational investigations into other dynamically important motions in DNA.

Atomic Mechanisms of Transthyretin Tetramer Dissociation Studied by Molecular Dynamics Simulations

The dissociation of the transthyretin (TTR) tetramer into a monomer is closely related to various TTR amyloidoses in humans. While the tetramer dissociation has been reported to be the rate-limiting step for TTR aggregation, few details are known about the mechanism. Here, molecular dynamics (MD) simulations were performed by combining conventional MD and biased metadynamics to investigate the mechanism for the wild-type (WT) and mutant (T119M) structures. Both were found to have a great deal in common. Conventional MD simulations reveal that interfacial hydrophobic interactions contribute significantly to stabilize the tetramer. Interfacial residues including L17, V20, L110, and V121 with close contacts form a hydrophobic channel. Metadynamics simulations indicate that the mouth opening of the hydrophobic channel is the first and the most difficult step for dissociation. Interactions of V20 between opposing dimers lock four monomers into the tetramer, and disruption of the interactions is found to be involved in the final step. During the dissociation, an increasing extent of solvation was observed by calculating the radial distribution functions of water around interfacial hydrophobic residues, suggesting that water plays a role in driving the tetramer dissociation. Moreover, compared to T119, residue M119 has a longer side chain that extends into the hydrophobic channel, making solvation more difficult, consistent with a higher energy barrier for dissociation of the T119M tetramer. This result provides a good explanation for the protective role of the T119M mutation. Overall, this study can provide atomic-level insights to better understand the pathogenesis of TTR amyloidosis and guide rational drug design in the future.

DeepCV: A Deep Learning Framework for Blind Search of Collective Variables in Expanded Configurational Space

We present Deep learning for Collective Variables (DeepCV), a computer code that provides an efficient and customizable implementation of the deep autoencoder neural network (DAENN) algorithm that has been developed in our group for computing collective variables (CVs) and can be used with enhanced sampling methods to reconstruct free energy surfaces of chemical reactions. DeepCV can be used to conveniently calculate molecular features, train models, generate CVs, validate rare events from sampling, and analyze a trajectory for chemical reactions of interest. We use DeepCV in an example study of the conformational transition of cyclohexene, where metadynamics simulations are performed using DAENN-generated CVs. The results show that the adopted CVs give free energies in line with those obtained by previously developed CVs and experimental results. DeepCV is open-source software written in Python/C++ object-oriented languages, based on the TensorFlow framework and distributed free of charge for noncommercial purposes, which can be incorporated into general molecular dynamics software. DeepCV also comes with several additional tools, i.e., an application program interface (API), documentation, and tutorials.

Unraveling the Abnormal Molecular Mechanism of Suicide Inhibition of Cytochrome P450 3A4

Suicide inhibition of the CYP3A4 enzyme by a drug inactivates the enzyme in the drug biotransformation process and often shows safety concerns about the drug. Despite extensive experimental studies, the abnormal molecular mechanism of a suicide inhibitor that forms a covalent bond with the residue far away from the catalytically active center of CYP3A4 inactivating the enzyme remains elusive. Here, the authors used molecular simulation approaches to study in detail how diquinone methide (DQR), the metabolite product of raloxifene, unbinds from CYP3A4 and inactivates the enzyme at the atomistic level. The results clearly indicate that in one of the intermediate states formed in its unbinding process, DQR covalently binds to Cys239, a residue far away from the catalytically active center of CYP3A4, and hinders the substrate from entering or leaving the enzyme. This work therefore provides an unprecedented way of clarifying the abnormal mechanism of suicide inhibition of the CYP3A4 enzyme.

Membrane Lipids Are an Integral Part of Transmembrane Allosteric Sites in GPCRs: A Case Study of Cannabinoid CB1 Receptor Bound to a Negative Allosteric Modulator, ORG27569, and Analogs

A growing number of G-protein-coupled receptor (GPCR) structures reveal novel transmembrane lipid-exposed allosteric sites. Ligands must first partition into the surrounding membrane and take lipid paths to these sites. Remarkably, a significant part of the bound ligands appears exposed to the membrane lipids. The experimental structures do not usually account for the surrounding lipids, and their apparent contribution to ligand access and binding is often overlooked and poorly understood. Using classical and enhanced molecular dynamics simulations, we show that membrane lipids are critical in the access and binding of ORG27569 and its analogs at the transmembrane site of cannabinoid CB1 receptor. The observed differences in the binding affinity and cooperativity arise from the functional groups that interact primarily with lipids. Our results demonstrate the significance of incorporating membrane lipids as an integral component of transmembrane sites for accurate characterization, binding-affinity calculations, and lead optimization in drug discovery.

Multiscale Simulations Identify Origins of Differential Carbapenem Hydrolysis by the OXA-48 β-Lactamase

OXA-48 β-lactamases are frequently encountered in bacterial infections caused by carbapenem-resistant Gram-negative bacteria. Due to the importance of carbapenems in the treatment of healthcare-associated infections and the increasingly wide dissemination of OXA-48-like enzymes on plasmids, these β-lactamases are of high clinical significance. Notably, OXA-48 hydrolyzes imipenem more efficiently than other commonly used carbapenems, such as meropenem. Here, we use extensive multiscale simulations of imipenem and meropenem hydrolysis by OXA-48 to dissect the dynamics and to explore differences in the reactivity of the possible conformational substates of the respective acylenzymes. Quantum mechanics/molecular mechanics (QM/MM) simulations of the deacylation reaction for both substrates demonstrate that deacylation is favored when the 6α-hydroxyethyl group is able to hydrogen bond to the water molecule responsible for deacylation but disfavored by the increasing hydration of either oxygen of the carboxylated Lys73 general base. Differences in free energy barriers calculated from the QM/MM simulations correlate well with the experimentally observed differences in hydrolytic efficiency between meropenem and imipenem. We conclude that the impaired breakdown of meropenem, compared to imipenem, which arises from a subtle change in the hydrogen bonding pattern between the deacylating water molecule and the antibiotic, is most likely induced by the meropenem 1β-methyl group. In addition to increased insights into carbapenem breakdown by OXA β-lactamases, which may aid in future efforts to design antibiotics or inhibitors, our approach exemplifies the combined use of atomistic simulations in determining the possible different enzyme-substrate substates and their influence on enzyme reaction kinetics.

Interconversion between Serum Amyloid A Native and Fibril Conformations

Overexpression of serum amyloid A (SAA) can lead to a form of amyloidosis where the fibrils are made of SAA fragments, most often SAA1-76. Using Replica Exchange with Tunneling, we study the conversion of a SAA1-76 chain between the folded conformation and a fibril conformation. We find that the basins in the free energy landscape corresponding to the two motifs are separated by barriers of only about 2-3 k B T. Crucial for the assembly into the fibril structure is the salt bridge 26E-34K that provides a scaffold for forming the fibril conformation.

Residue Folding Degree-Relationship to Secondary Structure Categories and Use as Collective Variable

Recently, we have shown that the residue folding degree, a network-based measure of folded content in proteins, is able to capture backbone conformational transitions related to the formation of secondary structures in molecular dynamics (MD) simulations. In this work, we focus primarily on developing a collective variable (CV) for MD based on this residue-bound parameter to be able to trace the evolution of secondary structure in segments of the protein. We show that this CV can do just that and that the related energy profiles (potentials of mean force, PMF) and transition barriers are comparable to those found by others for particular events in the folding process of the model mini protein Trp-cage. Hence, we conclude that the relative segment folding degree (the newly proposed CV) is a computationally viable option to gain insight into the formation of secondary structures in protein dynamics. We also show that this CV can be directly used as a measure of the amount of α-helical content in a selected segment.

Closer Look at Inverse Electron Demand Diels–Alder and Nucleophilic Addition Reactions on s-Tetrazines Using Enhanced Sampling Methods

Inverse electron demand [4+2] Diels–Alder (iEDDA) reactions as well as unprecedented nucleophilic (azaphilic) additions of R-substituted silyl-enol ethers (where R is Phenyl, Methyl, and Hydrogen) to 1,2,4,5-tetrazine (s-tetrazine) catalyzed by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {BF}_{3}$$\end{document}BF3 have recently been discovered (Simon et al. in Org Lett 23(7):2426–2430, 2021), where static calculations were employed for calculation of activation energies. In order to have a more realistic dynamic description of these reactions in explicit solution at ambient conditions, in this work we use a semiempirical tight-binding method combined with enhanced sampling techniques to calculate free energy surfaces (FESs) of the iEDDA and azaphilic addition reactions. Relevant products of not only s-tetrazine but also its derivatives such as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {BF}_{3}$$\end{document}BF3-mediated s-tetrazine adducts are investigated. We reconstruct the FESs of the iEDDA and azaphilic addition reactions using metadynamics and blue moon ensemble, and compare the ability of different collective variables (CVs) including bond distances, Social PeRmutation INvarianT (SPRINT) coordinates, and path-CV to describe the reaction pathway. We find that when a bulky Phenyl is used as a substituent at the dienophile the azaphilic addition is preferred over the iEDDA reaction. In addition, we also investigate the effect of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {BF}_{3}$$\end{document}BF3 in the diene and steric hindrance in the dienophile on the competition between the iEDDA and azaphilic addition reactions, providing chemical insight for reaction design.

Out-of-Equilibrium Biophysical Chemistry: The Case for Multidimensional, Integrated Single-Molecule Approaches

Out-of-equilibrium processes are ubiquitous across living organisms and all structural hierarchies of life. At the molecular scale, out-of-equilibrium processes (for example, enzyme catalysis, gene regulation, and motor protein functions) cause biological macromolecules to sample an ensemble of conformations over a wide range of time scales. Quantifying and conceptualizing the structure–dynamics to function relationship is challenging because continuously evolving multidimensional energy landscapes are necessary to describe nonequilibrium biological processes in biological macromolecules. In this perspective, we explore the challenges associated with state-of-the-art experimental techniques to understanding biological macromolecular function. We argue that it is time to revisit how we probe and model functional out-of-equilibrium biomolecular dynamics. We suggest that developing integrated single-molecule multiparametric force–fluorescence instruments and using advanced molecular dynamics simulations to study out-of-equilibrium biomolecules will provide a path towards understanding the principles of and mechanisms behind the structure–dynamics to function paradigm in biological macromolecules.

Independent Nontargeted Parallel Cascade Selection Molecular Dynamics (Ino-PaCS-MD) to Enhance the Conformational Sampling of Proteins

Biological functions are related to long-time protein dynamics (rare events) that are induced over microseconds. Such protein dynamics can be investigated using molecular dynamics (MD) simulations. However, the detection of rare events remains challenging using conventional MD (cMD) since the accessible timescales of cMD are shorter than those of the biological functions. Recently, the parallel cascade selection MD (PaCS-MD) has been proposed to detect such rare events, wherein transition paths are generated between a given reactant and product. As an extension, the nontargeted PaCS-MD (nt-PaCS-MD) has been proposed to predict the transition paths without requiring reference to any product. Thus, as a further extension, we herein propose independent nt-PaCS-MD, namely, Ino-PaCS-MD, wherein multiple walkers are launched from a set of different starting configurations. Each walker repeats a cycle of restarting short-time MD simulations from configurations with high potentials for making transitions to neighboring metastable states. To further enhance the sampling ability, Ino-PaCS-MD temporarily stops the conformational search and periodically resets the starting configurations so that they are uniformly distributed in a conformational subspace, thereby preventing a given protein from being trapped in one of the metastable states. As a demonstration, Ino-PaCS-MD successfully detects rare events of a maltose-binding protein as open-close transitions with a nanosecond-order simulation time, although a microsecond-order cMD simulation failed to detect these rare events, showing the high sampling efficiency of Ino-PaCS-MD.

Rate of Insulin Dimer Dissociation: Interplay between Memory Effects and Higher Dimensionality

We calculate the rate of dissociation of an insulin dimer into two monomers in water. The rate of this complex reaction is determined by multiple factors that are elucidated. By employing advanced sampling techniques, we first obtain the reaction free energy surface for the dimer dissociation as a function of two order parameters, namely, the distance between the center-of-mass of two monomers (R) and the number of cross-contacts (Q) among the backbone C_α atoms of two monomers. We then construct an orthogonal 2D reaction energy surface by introducing the reaction coordinate X to denote the minimum energy pathway and a conjugate coordinate Y that spans the orthogonal direction. The free energy landscape is rugged with multiple maxima and minima. We calculate the rate by employing not only the non-Markovian multidimensional rate theory but also several other theoretical approaches. The necessary reaction frequencies and the frictions are calculated from the time correlation function formalism. Our best estimate of the rate is 0.4 μs^-1. Our study reveals interesting opposite influences of dimensionality and memory in determining the rate constant of the reaction. We gain interesting insights into the dimer dissociation process by looking directly at the trajectories obtained from molecular dynamics simulation.

Membrane-facilitated receptor access and binding mechanisms of long-acting β2-adrenergic receptor (β2-AR) agonists

The drugs salmeterol, formoterol, and salbutamol constitute the frontline treatment for asthma and other chronic pulmonary diseases. These drugs activate the β2-adrenergic receptors (β2-AR), a class A G-protein-coupled receptor (GPCR) and differ significantly in their clinical onset and duration of actions. According to the "microkinetic model," the long duration of action of salmeterol and formoterol compared to salbutamol were attributed, at least in part, to their high lipophilicity and increased local concentrations in the membrane near the receptor. However, the structural and molecular bases of how the lipophilic drugs reach the binding site of the receptor from the surrounding membrane remain unknown. Using a variety of classical and enhanced molecular dynamics simulation techniques, we investigated the membrane partitioning characteristics, binding, and unbinding mechanisms of the ligands. The obtained results offer remarkable insight into the functional role of membrane lipids in the ligand association process. Strikingly, salmeterol entered the binding site from the bilayer through transmembrane helices 1 and 7. The entry was preceded by membrane-facilitated rearrangement and presentation of its phenyl-alkoxy-alkyl tail as a passkey to an access route gated by F193, a residue known critical for salmeterol's affinity. Formoterol's access is through the aqueous path shared by other β2-AR agents. We observed a novel secondary path for salbutamol that is distinct from its primary route. Our study offers a mechanistic description for the membrane-facilitated access and binding of ligands to β2-AR and establishes a groundwork for recognizing membrane lipids as an integral component in the molecular recognition process. Significance Statement The cell membrane's functional role behind the duration of action of long-acting β2-adrenergic receptor (β2-AR) agonists such as salmeterol has been a subject of debate for a long time. We investigated the binding and unbinding mechanisms of the three commonly used β2-AR agonists, salmeterol, formoterol, and salbutamol, using advanced simulation techniques. The obtained results offer unprecedented insights into the active role of membrane lipids in facilitating access and binding of the ligands, affecting the molecular recognition process and their pharmacology.

Molecular Routes of Dynamic Autocatalysis for Methanol-to-Hydrocarbons Reaction

The industrially important methanol-to-hydrocarbons (MTH) reaction is driven and sustained by autocatalysis in a dynamic and complex manner. Hitherto, the entire molecular routes and chemical nature of the autocatalytic network have not been well understood. Herein, with a multitechnique approach and multiscale analysis, we have obtained a full theoretical picture of the domino cascade of autocatalytic reaction network taking place on HZSM-5 zeolite. The autocatalytic reaction is demonstrated to be plausibly initiated by reacting dimethyl ether (DME) with the surface methoxy species (SMS) to generate the initial olefins, as evidenced by combining the kinetic analysis, in situ DRIFT spectroscopy, 2D ¹³C-¹³C MAS NMR, electronic states, and projected density of state (PDOS) analysis. This process is operando tracked and visualized at the picosecond time scale by advanced ab initio molecular dynamics (AIMD) simulations. The initial olefins ignite autocatalysis by building the first autocatalytic cycle-olefins-based cycle-followed by the speciation of methylcyclopentenyl (MCP) and aromatic cyclic active species. In doing so, the active sites accomplish the dynamic evolution from proton acid sites to supramolecular active centers that are experimentally identified with an ever-evolving and fluid feature. The olefins-guided and cyclic-species-guided catalytic cycles are interdependently linked to forge a previously unidentified hypercycle, being composed of one "selfish" autocatalytic cycle (i.e., olefins-based cycle with lighter olefins as autocatalysts for catalyzing the formation of olefins) and three cross-catalysis cycles (with olefinic, MCP, and aromatic species as autocatalysts for catalyzing each other's formation). The unraveled dynamic autocatalytic cycles/network would facilitate the catalyst design and process control for MTH technology.

Unfolding of Dynamical Events in the Early Stage of Insulin Dimer Dissociation

The dissociation of an insulin dimer is an important biochemical event that could also serve as a prototype of dissociations in similar biomolecular assemblies. We use a recently developed multidimensional free energy landscape for insulin dimer dissociation to unearth the microscopic and mechanistic aspects of the initial stages of the process that could hold the key to understanding the stability and the rate. The following sequence of events occurs in the initial stages: (i) The backbone hydrogen bonds break partially at the antiparallel β-sheet junction, (ii) the two α-helices (chain B) move away from each other while several residues (chain A) move closer, and (iii) a flow of adjacent water molecules occurs into the junction region. Interestingly, the intermonomeric center-to-center distance does not increase, but the number of native contacts exhibits a sharp decrease. Subsequent steps involve further disengagement of hydrophobic groups. This process is slow because of an entropic bottleneck created by the existence of the large configuration space available in the native state (NS), which is inhabited by low-frequency conformational fluctuations. We carry out a density-of-states analyses in the dimer NS to unearth distinctive features not present in the monomers. These low-frequency modes are also responsible for a large entropic stabilization of the NS. Hydrophobic disengagement in the early stage leads to the formation of a twisted intermediate state which itself is a metastable minimum (IS-1). The subsequent progress leads to another dimeric complex (IS-2), which is on the dissociative pathway and characterized by a further decrease in the native contacts. The dissociation process provides insights into the workings of a biomolecular assembly.

Loop Dynamics and Enzyme Catalysis in Protein Tyrosine Phosphatases

Protein tyrosine phosphatases (PTPs) play an important role in cellular signaling and have been implicated in human cancers, diabetes, and obesity. Despite shared catalytic mechanisms and transition states for the chemical steps of catalysis, catalytic rates within the PTP family vary over several orders of magnitude. These rate differences have been implied to arise from differing conformational dynamics of the closure of a protein loop, the WPD-loop, which carries a catalytically critical residue. The present work reports computational studies of the human protein tyrosine phosphatase 1B (PTP1B) and YopH from Yersinia pestis, for which NMR has demonstrated a link between their respective rates of WPD-loop motion and catalysis rates, which differ by an order of magnitude. We have performed detailed structural analysis, both conventional and enhanced sampling simulations of their loop dynamics, as well as empirical valence bond simulations of the chemical step of catalysis. These analyses revealed the key residues and structural features responsible for these differences, as well as the residues and pathways that facilitate allosteric communication in these enzymes. Curiously, our wild-type YopH simulations also identify a catalytically incompetent hyper-open conformation of its WPD-loop, sampled as a rare event, previously only experimentally observed in YopH-based chimeras. The effect of differences within the WPD-loop and its neighboring loops on the modulation of loop dynamics, as revealed in this work, may provide a facile means for the family of PTP enzymes to respond to environmental changes and regulate their catalytic activities.

Two distinct catalytic pathways for GH43 xylanolytic enzymes unveiled by X-ray and QM/MM simulations

Xylanolytic enzymes from glycoside hydrolase family 43 (GH43) are involved in the breakdown of hemicellulose, the second most abundant carbohydrate in plants. Here, we kinetically and mechanistically describe the non-reducing-end xylose-releasing exo-oligoxylanase activity and report the crystal structure of a native GH43 Michaelis complex with its substrate prior to hydrolysis. Two distinct calcium-stabilized conformations of the active site xylosyl unit are found, suggesting two alternative catalytic routes. These results are confirmed by QM/MM simulations that unveil the complete hydrolysis mechanism and identify two possible reaction pathways, involving different transition state conformations for the cleavage of xylooligosaccharides. Such catalytic conformational promiscuity in glycosidases is related to the open architecture of the active site and thus might be extended to other exo-acting enzymes. These findings expand the current general model of catalytic mechanism of glycosidases, a main reaction in nature, and impact on our understanding about their interaction with substrates and inhibitors.

Increasing the Sampling Efficiency of Protein Conformational Change by Combining a Modified Replica Exchange Molecular Dynamics and Normal Mode Analysis

Understanding conformational change at an atomic level is significant when determining a protein functional mechanism. Replica exchange molecular dynamics (REMD) is a widely used enhanced sampling method to explore protein conformational space. However, REMD with an explicit solvent model requires huge computational resources, immensely limiting its application. In this study, a variation of parallel tempering metadynamics (PTMetaD) with the omission of solvent-solvent interactions in exchange attempts and the use of low-frequency modes calculated by normal-mode analysis (NMA) as collective variables (CVs), namely ossPTMetaD, is proposed with the aim to accelerate MD simulations simultaneously in temperature and geometrical spaces. For testing the performance of ossPTMetaD, five protein systems with diverse biological functions and motion patterns were selected, including large-scale domain motion (AdK), flap movement (HIV-1 protease and BACE1), and DFG-motif flip in kinases (p38α and c-Abl). The simulation results showed that ossPTMetaD requires much fewer numbers of replicas than temperature REMD (T-REMD) with a reduction of ∼70% to achieve a similar exchange ratio. Although it does not obey the detailed balance condition, ossPTMetaD provides consistent results with T-REMD and experimental data. The high accessibility of the large conformational change of protein systems by ossPTMetaD, especially in simulating the very challenging DFG-motif flip of protein kinases, demonstrated its high efficiency and robustness in the characterization of the large-scale protein conformational change pathway and associated free energy profile.

Bifurcated Hydrogen Bonds and the Fold Switching of Lymphotactin

Lymphotactin (Ltn) exists under physiological conditions in an equilibrium between two interconverting structures with distinct biological functions. Using replica-exchange-with-tunneling, we study the conversion between the 2-folds. Unlike previously proposed, we find that the fold switching does not require unfolding of lymphotactin but proceeds through a series of intermediates that remain partially structured. This process relies on two bifurcated hydrogen bonds that connect the β2 and β3 strands and ease the transition between the hydrogen bond pattern by which the central three-stranded β-sheet in the two forms differs.

POLARIS: Path of Least Action Analysis on Energy Landscapes

Free-energy landscapes are a powerful tool for analyzing dynamical processes - capable of providing a complete mapping of a system's configurations in state space while articulating its energetics topologically in the form of sprawling hills and valleys. Within this mapping, the path of least action can be derived - representing the most probable sequence of transitions taken between any two states in the landscape. In this article, POLARIS (Path of Least Action Recursive Survey) is presented as a dynamic, global approach that efficiently automates the discovery of the least action path on previously determined 2D energy landscapes. Important built-in features of this program include plotting of landscape trajectories and transition state theory diagrams, generation of text files with least action coordinates and respective energies, and bifurcation analysis tools that provide downstream versatility for comparing most probable paths and reaction rates.

Drivers of α-Sheet Formation in Transthyretin under Amyloidogenic Conditions

Amyloid diseases make up a set of fatal disorders in which proteins aggregate to form fibrils that deposit in tissues throughout the body. Amyloid-associated diseases are challenging to study because amyloid formation occurs on time scales that span several orders of magnitude and involve heterogeneous, interconverting protein conformations. The development of more effective technologies to diagnose and treat amyloid disease requires both a map of the conformations sampled during amyloidogenesis and an understanding of the molecular mechanisms that drive this process. In prior molecular dynamics simulations of amyloid proteins, we observed the formation of a nonstandard type of secondary structure, called α-sheet, that we proposed is associated with the pathogenic conformers in amyloid disease, the soluble oligomers. However, the detailed molecular interactions that drive the conversion to α-sheet remain elusive. Here we use molecular dynamics simulations to interrogate a critical event in transthyretin aggregation, the formation of aggregation-competent, monomeric species. We show that conformational changes in one of the two β-sheets in transthyretin enable solvent molecules and polar side chains to form electrostatic interactions with main-chain peptide groups to facilitate and modulate conversion to α-sheet secondary structure. Our results shed light on the early conformational changes that drive transthyretin toward the α-sheet structure associated with toxicity. Delineation of the molecular events that lead to aggregation at atomic resolution can aid strategies to target the early, critical toxic soluble oligomers.

Exploring the Ligand Binding/Unbinding Pathway by Selectively Enhanced Sampling of Ligand in a Protein-Ligand Complex

Understanding the protein-ligand binding is of fundamental biological interest and is essential for structure-based drug design. The difficulty in capturing the dynamic process, however, poses a great challenge for current experimental and theoretical simulation techniques. A selective integrated-tempering-sampling molecular dynamics (SITSMD) method offering an option for selectively enhanced sampling of the ligand in a protein-ligand complex was utilized to quantitatively illuminate the binding of benzamidine to the wild-type trypsin protease and its two mutants (S214E and S214K). The SITSMD simulations could produce consistent results as the extensive conventional MD simulation and gave additional insights into the binding pathway for the test protein-ligand complex system using significantly saved computational resource and time, indicating the potential of such a method in investigating protein-ligand binding. Additionally, the simulations identified the different roles of trypsin-benzamidine van der Waals (vdW) and electrostatic interactions in the binding: the former interaction works as the driving force for dragging the benzamidine close to the native binding pocket, and the latter interaction mainly contributes to stabilizing the benzamidine inside the pocket. The S214E mutation introduces more favorable electrostatic interactions, and as a result, both vdW and electrostatic interactions drive the benzamidine binding, lowering the binding and unbinding free energy barrier. In contrast, the S214K mutation prohibits the binding of the benzamidine to the native ligand binding pocket by introducing disliked charge-charge interactions. In summary, these findings suggest that the change in specific residues could modify the protein druggability, including the binding kinetics and thermodynamics.