All-atom empirical potential for molecular modeling and dynamics studies of proteins

Identification and Ranking of Binding Sites from Structural Ensembles: Application to SARS-CoV-2

Target identification and evaluation is a critical step in the drug discovery process. Although time-intensive and complex, the challenge becomes even more acute in the realm of infectious disease, where the rapid emergence of new viruses, the swift mutation of existing targets, and partial effectiveness of approved antivirals can lead to outbreaks of significant public health concern. The COVID-19 pandemic, caused by the SARS-CoV-2 virus, serves as a prime example of this, where despite the allocation of substantial resources, Paxlovid is currently the only effective treatment. In that case, significant effort pre-pandemic had been expended to evaluate the biological target for the closely related SARS-CoV. In this work, we utilize the computational hot spot mapping method, FTMove, to rapidly identify and rank binding sites for a set of nine SARS-CoV-2 drug/potential drug targets. FTMove takes into account protein flexibility by mapping binding site hot spots across an ensemble of structures for a given target. To assess the applicability of the FTMove approach to a wide range of drug targets for viral pathogens, we also carry out a comprehensive review of the known SARS-CoV-2 ligandable sites. The approach is able to identify the vast majority of all known sites and a few additional sites, which may in fact be yet to be discovered as ligandable. Furthermore, a UMAP analysis of the FTMove features for each identified binding site is largely able to separate predicted sites with experimentally known binders from those without known binders. These results demonstrate the utility of FTMove to rapidly identify actionable sites across a range of targets for a given indication. As such, the approach is expected to be particularly useful for assessing target binding sites for any emerging pathogen, as well as for indications in other disease areas, and providing actionable starting points for structure-based drug design efforts.

A Deep Learning-Driven Sampling Technique to Explore the Phase Space of an RNA Stem-Loop

The folding and unfolding of RNA stem-loops are critical biological processes; however, their computational studies are often hampered by the ruggedness of their folding landscape, necessitating long simulation times at the atomistic scale. Here, we adapted DeepDriveMD (DDMD), an advanced deep learning-driven sampling technique originally developed for protein folding, to address the challenges of RNA stem-loop folding. Although tempering- and order parameter-based techniques are commonly used for similar rare-event problems, the computational costs or the need for a priori knowledge about the system often present a challenge in their effective use. DDMD overcomes these challenges by adaptively learning from an ensemble of running MD simulations using generic contact maps as the raw input. DeepDriveMD enables on-the-fly learning of a low-dimensional latent representation and guides the simulation toward the undersampled regions while optimizing the resources to explore the relevant parts of the phase space. We showed that DDMD estimates the free energy landscape of the RNA stem-loop reasonably well at room temperature. Our simulation framework runs at a constant temperature without external biasing potential, hence preserving the information on transition rates, with a computational cost much lower than that of the simulations performed with external biasing potentials. We also introduced a reweighting strategy for obtaining unbiased free energy surfaces and presented a qualitative analysis of the latent space. This analysis showed that the latent space captures the relevant slow degrees of freedom for the RNA folding problem of interest. Finally, throughout the manuscript, we outlined how different parameters are selected and optimized to adapt DDMD for this system. We believe this compendium of decision-making processes will help new users adapt this technique for the rare-event sampling problems of their interest.

Potential toxicity of Graphene (Oxide) quantum dots to human intestinal fatty acid binding protein (HIFABP) via obstructing the protein's openings

Graphene quantum dots (GQDs) have garnered significant attention across numerous fields due to their ultrasmall size and exceptional properties. However, their extensive applications may lead to environmental exposure and subsequent uptake by humans. Yet, conflicting reports exist regarding the potential toxicity of GQDs based on experimental investigations. Therefore, a comprehensive understanding of GQD biosafety requires further microscopic and molecular-level investigations. In this study, we employed molecular dynamics (MD) simulations to explore the interactions between GQDs and graphene oxide quantum dots (GOQDs) with a protein model, the human intestinal fatty acid binding protein (HIFABP), that plays a crucial role in mediating the carrier of fatty acids in the intestine. Our MD simulation results reveal that GQDs can be adsorbed on the opening of HIFABP, which serves as an entrance for the fatty acid molecules into the protein's interior cavity. This adsorption has the potential to obstruct the opening of HIFABP, leading to the loss of its normal biological function and ultimately resulting in toxicity. The adsorption of GQDs is driven by a combination of van der Waals (vdW), π-π stacking, cation-π, and hydrophobic interactions. Similarly, GOQDs also exhibit the ability to block the opening of HIFABP, thereby potentially causing toxicity. The blockage of GOQDs to HIFABP is guided by a combination of vdW, Coulomb, π-π stacking, and hydrophobic interactions. These findings not only highlight the potential harmful effects of GQDs on HIFABP but also elucidate the underlying molecular mechanism, which provides crucial insights into GQD toxicology.

The Generation of ROS by Exposure to Trihalomethanes Promotes the IκBα/NF-κB/p65 Complex Dissociation in Human Lung Fibroblast

Background: Disinfection by-products used to obtain drinking water, including halomethanes (HMs) such as CH2Cl2, CHCl3, and BrCHCl2, induce cytotoxicity and hyperproliferation in human lung fibroblasts (MRC-5). Enzymes such as superoxide dismutase (SOD), catalase (CAT), and glutathione peroxidase (GPx) modulate these damages through their biotransformation processes, potentially generating toxic metabolites. However, the role of the oxidative stress response in cellular hyperproliferation, modulated by nuclear factor-kappa B (NF-κB), remains unclear. Methods: In this study, MRC-5 cells were treated with these compounds to evaluate reactive oxygen species (ROS) production, lipid peroxidation, phospho-NF-κB/p65 (Ser536) levels, and the activities of SOD, CAT, and GPx. Additionally, the interactions between HMs and ROS with the IκBα/NF-κB/p65 complex were analyzed using molecular docking. Results: Correlation analysis among biomarkers revealed positive relationships between pro-oxidant damage and antioxidant responses, particularly in cells treated with CH2Cl2 and BrCHCl2. Conversely, negative relationships were observed between ROS levels and NF-κB/p65 levels in cells treated with CH2Cl2 and CHCl3. The estimated relative free energy of binding using thermodynamic integration with the p65 subunit of NF-κB was -3.3 kcal/mol for BrCHCl2, -3.5 kcal/mol for both CHCl3 and O2•, and -3.6 kcal/mol for H2O2. Conclusions: Chloride and bromide atoms were found in close contact with IPT domain residues, particularly in the RHD region involved in DNA binding. Ser281 is located within this domain, facilitating the phosphorylation of this protein. Similarly, both ROS interacted with the IPT domain in the RHD region, with H2O2 forming a side-chain oxygen interaction with Leu280 adjacent to the phosphorylation site of p65. However, the negative correlation between ROS and phospho-NF-κB/p65 suggests that steric hindrance by ROS on the C-terminal domain of NF-κB/p65 may play a role in the antioxidant response.

Steviol rebaudiosides bind to four different sites of the human sweet taste receptor (T1R2/T1R3) complex explaining confusing experiments

Sucrose provides both sweetness and energy by binding to both Venus flytrap domains (VFD) of the heterodimeric sweet taste receptor (T1R2/T1R3). In contrast, non-caloric sweeteners such as sucralose and aspartame only bind to one specific domain (VFD2) of T1R2, resulting in high-intensity sweetness. In this study, we investigate the binding mechanism of various steviol glycosides, artificial sweeteners, and a negative allosteric modulator (lactisole) at four distinct binding sites: VFD2, VFD3, transmembrane domain 2 (TMD2), and TMD3 through binding experiments and computational docking studies. Our docking results reveal multiple binding sites for the tested ligands, including the radiolabeled ligands. Our experimental evidence demonstrates that the C20 carboxy terminus of the Gα protein can bind to the intracellular region of either TMD2 or TMD3, altering GPCR affinity to the high-affinity state for steviol glycosides. These findings provide a mechanistic understanding of the structure and function of this heterodimeric sweet taste receptor.

Dissecting the neuroprotective interaction between the BH4 domain of BCL-w and the IP3 receptor

BCL-w is a BCL-2 family protein that promotes cell survival in tissue- and disease-specific contexts. The canonical anti-apoptotic functionality of BCL-w is mediated by a surface groove that traps the BCL-2 homology 3 (BH3) α-helices of pro-apoptotic members, blocking cell death. A distinct N-terminal portion of BCL-w, termed the BCL-2 homology 4 (BH4) domain, selectively protects axons from paclitaxel-induced degeneration by modulating IP3 receptors, a noncanonical BCL-2 family target. Given the potential of BCL-w BH4 mimetics to prevent or mitigate chemotherapy-induced peripheral neuropathy, we sought to characterize the interaction between BCL-w BH4 and the IP3 receptor, combining "staple" and alanine scanning approaches with molecular dynamics simulations. We generated and identified stapled BCL-w BH4 peptides with optimized IP3 receptor binding and neuroprotective activities. Point mutagenesis further revealed the sequence determinants for BCL-w BH4 specificity, providing a blueprint for therapeutic targeting of IP3 receptors to achieve neuroprotection.

CHARMM at 45: Enhancements in Accessibility, Functionality, and Speed

Since its inception nearly a half century ago, CHARMM has been playing a central role in computational biochemistry and biophysics. Commensurate with the developments in experimental research and advances in computer hardware, the range of methods and applicability of CHARMM have also grown. This review summarizes major developments that occurred after 2009 when the last review of CHARMM was published. They include the following: new faster simulation engines, accessible user interfaces for convenient workflows, and a vast array of simulation and analysis methods that encompass quantum mechanical, atomistic, and coarse-grained levels, as well as extensive coverage of force fields. In addition to providing the current snapshot of the CHARMM development, this review may serve as a starting point for exploring relevant theories and computational methods for tackling contemporary and emerging problems in biomolecular systems. CHARMM is freely available for academic and nonprofit research at https://academiccharmm.org/program.

Physics-Based Protein Networks Might Recover Effectful Mutations─a Case Study on Cathepsin G

Molecular dynamics simulations have been remarkably effective for observing and analyzing structures and dynamics of proteins, with longer trajectories being computed every day. Still, often, relevant time scales are not observed. Adequately analyzing the generated trajectories can highlight the interesting areas within a protein such as mutation sites or allosteric hotspots, which might foreshadow dynamics untouched by the simulations. We employ a physics-based protein network and propose that such a network can adequately analyze the protein dynamics. The analysis is conducted on simulations of cathepsin G and neutrophil elastase, which are remarkably similar but with different specificities. However, a single mutation in cathepsin G recovers the specificity of neutrophil elastase. The physics-based network built on the interactions between residues instead of the distances can pinpoint the active triad in the proteins studied. Overall, the network seems to capture the structural behavior better than purely distance-based networks.

Statistical-Mechanics Analyses on Thermodynamics of Protein Folding Constructed by Privalov and Co-Workers

Privalov and co-workers estimated the changes in hydration enthalpy and entropy upon ubiquitin unfolding and their temperature dependences denoted by ΔHhyd(T) and ΔShyd(T), respectively, from experimentally measured enthalpies and entropies of transfer of various model compounds from gaseous phase to water. We calculate ΔHhyd(T) and ΔShyd(T) for ubiquitin by our statistical-mechanics theory where molecular and atomistic models are employed for water and protein structure, respectively. ΔHhyd(T) and ΔShyd(T) calculated are in remarkably good agreement with those estimated by Privalov and co-workers. By examining relative magnitudes and signs of the changes in a variety of constituents of ΔHhyd(T) and ΔShyd(T), we confirm that the hydrophobic effect is an essential force driving a protein to fold. Detailed and comprehensive explanations are given for our claim that the prevailing views of the hydrophobic effect are not capable of elucidating its weakening at low temperatures, whereas our updated view is. We find out problematic points of the changes in enthalpy and entropy upon protein unfolding denoted by ΔH°(T) and ΔS°(T), respectively, which are measured using the differential scanning calorimetry at low pH, suggesting a theoretical method of calculating ΔH°(T) and ΔS°(T) at pH ∼ 7.

Ionizable Cationic Lipids and Helper Lipids Synergistically Contribute to RNA Packing and Protection in Lipid-Based Nanomaterials

Lipid-based nanomaterials are used as a common delivery vehicle for RNA therapeutics. They typically include a formulation containing ionizable cationic lipids, cholesterol, phospholipids, and a small molar fraction of PEGylated lipids. The ionizable cationic lipids are considered a crucial element of the formulation for the way they mediate interactions with the anionic RNA as a function of pH. Here, we show, by means of molecular dynamics simulation of lipid formulations containing two different ionizable cationic lipids (DLinDMA and DLinDAP), that the direct interactions of those lipids with RNA, taken alone, may not be sufficient to determine the level of protection and packaging of mRNA. Our simulations help and highlight how the collective behavior of the lipids in the formulation, which determines the ability to envelop the RNA, and the level of hydration of the lipid-RNA interface may also play a significant role. This allows the drawing of a hypothesis about the experimentally observed differences in the transfection efficiency of the two ionizable cationic lipids.

Mechanistic insights into carbonic anhydrase IX inhibition by coumarins from Calendula officinalis: in vitro and in silico approaches

Given the critical role of carbonic anhydrase IX (CA IX) in various pathological conditions, there is a significant demand for new inhibitors to enhance patient outcomes and clinical management. In this study, we examined the inhibitory effectiveness of five coumarins derived from Calendula officinalis against CA IX using in vitro assays and computational modeling. Among the coumarins tested, xeroboside and isobaisseoside were identified as the most potent inhibitors. Kinetic studies indicated that xeroboside and isobaisseoside exhibit a mixed inhibition mode. Molecular docking analysis showed that the tested coumarins exhibit binding affinities and extensive polar interactions with CA IX. These coumarins demonstrated significant hydrophobic interactions and occupied the same binding site as acetazolamide (AAZ). Molecular dynamics (MD) indicated that xeroboside and isobaisseoside exhibited consistent trajectories and notable energy stabilization during their interaction with CA IX. MM/PBSA calculations showed that xeroboside displayed the lowest binding free energy (-27.26 ± 2.48 kJ mol-1). Potential Energy Landscape (PEL) analysis revealed distinct and stable conformational states for the CA IX-ligand complexes, with xeroboside exhibiting the most stable and lowest energy configuration. These computational findings are consistent with the experimental results, highlighting the potential efficacy of xeroboside and isobaisseoside as CA IX inhibitors. In conclusion, Calendula officinalis-derived coumarins are promising candidates as effective CA IX inhibitors.

Plastocyanin and Cytochrome f Complex Structures Obtained by NMR, Molecular Dynamics, and AlphaFold 3 Methods Compared to Cryo-EM Data

Plastocyanin is a small mobile protein that facilitates electron transfer through the formation of short-lived protein-protein complexes with cytochrome bf and photosystem 1. Due to the transient nature of plastocyanin-cytochrome f complex, the lack of a long-lived tight complex makes it impossible to determine its structure by X-ray diffraction analysis. Up to today, a number of slightly different structures of such complexes have been obtained by experimental and computer methods. Now, artificial intelligence gives us the possibility to predict the structures of intermolecular complexes. In this study, we compare encounter and final complexes obtained by Brownian and molecular dynamics methods, as well as the structures predicted by AlphaFold 3, with NMR and cryo-EM data. Surprisingly, the best match for the plastocyanin electron density obtained by cryo-EM was demonstrated by an AlphaFold 3 structure. The orientation of plastocyanin in this structure almost completely coincides with its orientation obtained by molecular dynamics calculation, and, at the same time, it is different from the orientation of plastocyanin predicted on the basis of NMR data. This is even more unexpected given that only NMR structures for the plastocyanin-cytochrome f complex are available in the PDB database, which was used to train AlphaFold 3.

Dynamics-based protein network features accurately discriminate neutral and rheostat positions

In some proteins, a unique class of nonconserved positions is characterized by their ability to generate diverse functional outcomes through single amino acid substitutions. Due to their ability to tune protein function, accurately identifying such "rheostat" positions is crucial for protein design, for understanding the impact of mutations observed in humans, and for predicting the evolution of pathogen drug resistance. However, identifying rheostat positions has been challenging, due-in part-to the absence of a clear structural relationship with binding sites. In this study, experimental data from our previous study of the Escherichia coli lactose repressor protein (LacI) was used to identify rheostat positions for which mutations tune in vivo EC50 for the allosteric ligand "IPTG." We next used the rheostat assignments to test the hypothesis that rheostat positions have unique dynamic features that will enable their identification. To that end, we integrated all-atom molecular dynamics simulations with perturbation residue response analysis. Results first revealed distinct dynamic behavior in IPTG-bound LacI compared with apo LacI, which was consistent with IPTG's role as an allosteric inducer. Next, we used a variety of dynamic features to build a classification model that discriminates experimentally characterized rheostat positions in LacI from positions with other types of substitution outcomes. In parallel, we built a second classifier model based on the 3D structural "static" network features of LacI. In comparative studies, the dynamic model better identified rheostat positions that were >8 Å from the binding site. In summary, our study provides insights into the dynamic characteristics of rheostat positions and suggests that models built on dynamic features may be useful for predicting the locations of rheostat positions in a wide range of proteins.

Adenosine Triphosphate Inhibits Cold-Responsive Aggregation

Adenosine triphosphate (ATP), ubiquitous in all living organisms, is conventionally recognized as a fundamental energy currency essential for a myriad of cellular processes. While its traditional role in energy metabolism requires only micromolar concentrations, the cellular content of ATP has been found to be significantly higher at the millimolar level. Recent studies have attempted to correlate this higher concentration of ATP with its nonenergetic role in maintaining protein homeostasis, leaving the investigation of ATP's nontrivial activities in biology an open question. Here, by coupling computer simulations and experiments, we uncover new insights into ATP's role as a cryoprotectant against cold-salt stress, highlighting the necessity for higher cellular ATP concentrations. We present direct evidence at charged silica interfaces, demonstrating ATP's ability to restore native intersurface interactions disrupted by combined cold-salt stress, thereby inhibiting cold-responsive aggregation in high-salt conditions. ATP desorbs salt cations from negatively charged surfaces through predominant interactions between ATP and the salt cations. Although the mode of ATP's action remains unchanged with temperature, the extent of interaction scales with temperature, requiring less ATP activity at lower temperatures, justifying the reason for reduction in cellular ATP content due to the cold effect, reported in previous experimental studies. The trend observed in inorganic nanostructures is recurrent and robustly transferable to charged protein interfaces. A thorough comparison of ATP's cryoprotective activity with traditionally known biological cryoprotectants (glycine and betaine) reveals ATP's greater efficiency. In retrospect, our findings highlight ATP's additional biological role in cryopreservation, expanding its potential biomedical applications by offering effective protection of cells from cryoinjuries and avoiding the significant challenges associated with the toxicity of organic cryoprotectants.

Small Molecule NS11021 Promotes BK Channel Activation by Increasing Inner Pore Hydration

The Ca2+ and voltage-gated big potassium (BK) channels are implicated in various diseases, including heart disease, asthma, epilepsy, and cancer, but remain an elusive drug target. A class of negatively charged activators (NCAs) have been demonstrated to promote the activation of several potassium channels including BK channels by binding to the hydrophobic inner pore, yet the underlying molecular mechanism of action remains poorly understood. In this work, we analyze the binding mode and potential activation mechanism of a specific NCA named NS11021 using atomistic simulations. The results show that NS11021 binding to the pore in deactivated BK channels is nonspecific and dynamic. The binding free energy of -8.3 ± 0.7 kcal/mol (KD = 0.3-3.1 μM) calculated using umbrella sampling agrees quantitatively with the experimental EC50 range of 0.4-2.1 μM. The bound NS11021 remains dynamic and is distal from the filter to significantly impact its conformation. Instead, NS11021 binding significantly enhances the pore hydration due to the charged tetrazol-phenyl group, thereby promoting the opening of the hydrophobic gate. We further show that the free energy barrier to K+ permeation is reduced by ∼3 kcal/mol regardless of the binding pose, which could explain the ∼62-fold increase in the intrinsic opening of BK channels measured experimentally. Taken together, these results support the idea that the molecular mechanism of NS11021 derives from increasing the hydration level of the conformationally closed pore, which does not depend on specific binding and likely explains the ability of NCAs to activate multiple K+ channels.

Theoretical Investigation into Polymorphic Transformation between β-HMX and δ-HMX by Finite Temperature String

Polymorphic transformation is important in chemical industries, in particular, in those involving explosive molecular crystals. However, due to simulating challenges in the rare event method and collective variables, understanding the transformation mechanism of molecular crystals with a complex structure at the molecular level is poor. In this work, with the constructed order parameters (OPs) and K-means clustering algorithm, the potential of mean force (PMF) along the minimum free-energy path connecting β-HMX and δ-HMX was calculated by the finite temperature string method in the collective variables (SMCV), the free-energy profile and nucleation kinetics were obtained by Markovian milestoning with Voronoi tessellations, and the temperature effect on nucleation was also clarified. The barriers of transformation were affected by the finite-size effects. The configuration with the lower potential barrier in the PMF corresponded to the critical nucleus. The time and free-energy barrier of the polymorphic transformation were reduced as the temperature increased, which was explained by the pre-exponential factor and nucleation rate. Thus, the polymorphic transformation of HMX could be controlled by the temperatures, as is consistent with previous experimental results. Finally, the HMX polymorph dependency of the impact sensitivity was discussed. This work provides an effective way to reveal the polymorphic transformation of the molecular crystal with a cyclic molecular structure, and further to prepare the desired explosive by controlling the transformation temperature.

Diffusion and Viscosity in Mixed Protein Solutions

The viscosity and diffusion properties of crowded protein systems were investigated with molecular dynamics simulations of SH3 mixtures with different crowders, and results were compared with experimental data. The simulations accurately reproduced experimental trends across a wide range of protein concentrations, including highly crowded environments up to 300 g/L. Notably, viscosity increased with crowding but varied little between different crowder types, while diffusion rates were significantly reduced depending on protein-protein interaction strength. Analysis using the Stokes-Einstein relation indicated that the reduction in diffusion exceeded what was expected from viscosity changes alone, with the additional slow-down attributable to transient cluster formation driven by weakly attractive interactions. Contact kinetics analysis further revealed that longer-lived interactions contributed more significantly to reduced diffusion rates than short-lived interactions. This study also highlights the accuracy of current computational methodologies for capturing the dynamics of proteins in highly concentrated solutions and provides insights into the molecular mechanisms affecting protein mobility in crowded environments.

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

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

Structural insights into inhibitor mechanisms on immature HIV-1 Gag lattice revealed by high-resolution in situ single-particle cryo-EM

HIV-1 inhibitors, such as Bevirimat (BVM) and Lenacapavir (LEN), block the production and maturation of infectious virions. However, their mechanisms remain unclear due to the absence of high-resolution structures for BVM complexes and LEN's structural data being limited to the mature capsid. Utilizing perforated virus-like particles (VLPs) produced from mammalian cells, we developed an approach to determine in situ cryo-electron microscopy (cryo-EM) structures of HIV-1 with inhibitors. This allowed for the first structural determination of the native immature HIV-1 particle with BVM and LEN bound inside the VLPs at high resolutions. Our findings offer a more accurate model of BVM engaging the Gag lattice and, importantly, demonstrate that LEN not only binds the mature capsid but also targets the immature lattice in a distinct manner. The binding of LEN induces a conformational change in the capsid protein (CA) region and alters the architecture of the Gag lattice, which may affect the maturation process. These insights expand our understanding of the inhibitory mechanisms of BVM and LEN on HIV-1 and provide valuable clues for the design of future inhibitors.

Enhanced Sampling of Biomolecular Slow Conformational Transitions Using Adaptive Sampling and Machine Learning

Biomolecular simulations often suffer from the "time scale problem", hindering the study of rare events occurring over extended time scales. Enhanced sampling techniques aim to alleviate this issue by accelerating conformational transitions, yet they typically necessitate well-defined collective variables (CVs), posing a significant challenge. Machine learning offers promising solutions but typically requires rich training data encompassing the entire free energy surface (FES). In this work, we introduce an automated iterative pipeline designed to mitigate these limitations. Our protocol first utilizes a CV-free count-based adaptive sampling method to generate a data set rich in rare events. From this data set, slow modes are identified using Koopman-reweighted time-lagged independent component analysis (KTICA), which are subsequently leveraged by on-the-fly probability enhanced sampling (OPES) to efficiently explore the FES. The effectiveness of our pipeline is demonstrated and further compared with the common Markov State Model (MSM) approach on two model systems with increasing complexity: alanine dipeptide (Ala2) and deca-alanine (Ala10), underscoring its applicability across diverse biomolecular simulations.

Buckets Instead of Umbrellas for Enhanced Sampling and Free Energy Calculations

Umbrella sampling has been a workhorse for free energy calculations in molecular simulations for several decades. In conventional umbrella sampling, restraining bias potentials are strategically applied along one or several collective variables. Major drawbacks associated with this method are the requirement of a large number of bias windows and the poor sampling of the transverse coordinates. In this work, we propose an alternate formalism that departs from the traditional umbrella sampling to mitigate these issues, where we replace umbrella-type restraining bias potentials with bucket-type wall potentials. This modification permits one to formulate an efficient computational strategy leveraging wall potentials and metadynamics sampling. This new method, called "bucket sampling", can significantly reduce the computational cost of obtaining converged high-dimensional free energy surfaces. Extensions of the proposed method with temperature acceleration and replica exchange solute tempering are also demonstrated.

Investigating the Mechanisms of Antibody Binding to Alpha-Synuclein for the Treatment of Parkinson's Disease

Parkinson's disease (PD) is an idiopathic neurodegenerative disorder with the second-highest prevalence rate behind Alzheimer's disease. The pathophysiological hallmarks of PD are both degeneration of dopaminergic neurons in the substantia nigra pars compacta and the inclusion of misfolded α-synuclein (α-syn) aggregates known as Lewy bodies. Despite decades of research for potential PD treatments, none have been developed, and developing new therapeutic agents is a time-consuming and expensive process. Computational methods can be used to investigate the properties of drug candidates currently undergoing clinical trials to determine their theoretical efficiency at targeting α-syn. Monoclonal antibodies (mAbs) are biological drugs with high specificity, and Prasinezumab (PRX002) is an mAb currently in Phase II, which targets the C-terminus (AA 118-126) of α-syn. We utilized BioLuminate and PyMol for the structure prediction and preparation of the fragment antigen-binding (Fab) region of PRX002 and 34 different conformations of α-syn. Protein-protein docking simulations were performed using PIPER, and 3 of the docking poses were selected based on the best fit. Molecular dynamics simulations were conducted on the docked protein structures in triplicate for 1000 ns, and hydrogen bonds and electrostatic and hydrophobic interactions were analyzed using MDAnalysis to determine which residues were interacting and how often. Hydrogen bonds were shown to form frequently between the HCDR2 region of PRX002 and α-syn. Free energy was calculated to determine the binding affinity. The predicted binding affinity shows a strong antibody-antigen attraction between PRX002 and α-syn. RMSD was calculated to determine the conformational change of these regions throughout the simulation. The mAb's developability was determined using computational screening methods. Our results demonstrate the efficiency and developability of this therapeutic agent.

Revealing the interaction mechanism between bovine serum albumin (BSA) and a fluorescent coumarin derivative: A multispectroscopic and in silico approach

The interaction of biomacromolecules with each other or with the ligands is essential for biological activity. In this context, the molecular recognition of bovine serum albumin (BSA) with 4-(Benzo[1,3]dioxol-5-yloxymethyl)-7-hydroxy-chromen-2-one (4BHC) is explored using multispectroscopic and computational techniques. UV-Vis spectroscopy helped in predicting the conformational variations in BSA. Using fluorescence spectroscopy, the quenching behaviour of the fluorophore upon interaction with the ligand is examined, which is found to be a static type of quenching; fluorescence lifetime studies further verify this. The binding constant is discovered to be in the range of 104 M-1, which indicates the moderate type of association that results in reversible binding, where the transport and release of ligands in the target tissue takes place. Fourier-transform infrared spectroscopy (FT-IR) measurements validate the secondary structure conformational changes of BSA after complexing with 4BHC. The thermodynamic factors obtained through temperature-dependent fluorescence studies suggest that the dominant kind of interaction force is hydrophobic in nature, and the interaction process is spontaneous. The alterations in the surrounding microenvironment of the binding site and conformational shifts in the structure of the protein are studied through 3D fluorescence and synchronous fluorescence studies. Molecular docking and molecular dynamics (MD) simulations agree with experimental results and explain the structural stability throughout the discussion. The outcomes might have possible applications in the field of pharmacodynamics and pharmacokinetics.

Introductory Tutorials for Simulating Protein Dynamics with GROMACS

Atomistic molecular dynamics (MD) simulations have become an indispensable tool for investigating the structure, dynamics, and energetics of biomolecules. Continual optimization of software algorithms and hardware has enabled investigators to access biologically relevant time scales in feasible amounts of computing time. Given the widespread use and utility of MD simulations, there is considerable interest in learning essential skills in performing them. Here, we present a set of introductory tutorials for performing MD simulations of proteins in the popular, open-source GROMACS package. Three exercises are detailed, including simulating a single protein, setting up a protein complex, and performing umbrella sampling simulations to model the unfolding of a short polypeptide. Essential features and input settings are illustrated throughout. The purpose of these tutorials is to provide new users with a general understanding of foundational workflows, from which they can design their own simulations.

Osmolyte-Induced Modulation of Hofmeister Series

Natural selection has driven the convergence toward a selected set of osmolytes, endowing them with the necessary efficiency to manage stress arising from salt diversity. This study combines atomistic simulations and experiments to investigate how two osmolytes, glycine and betaine, individually modulate the Hofmeister ion ordering of alkali metal salts (LiCl, KCl, and CsCl) near a charged silica interface. Both osmolytes are found to prevent salt-induced aggregation of the charged entities, yet their mode and degree of relative modulation depend on their intricate interplay with specific salt cations. Betaine's ion-mediated surface interaction maintains Hofmeister ion ordering, whereas glycine alters the relative Hofmeister order of the cation by salt-specific ion desorption from the surface. Experimental validation through surface-enhanced Raman spectroscopy supports these findings, elucidating osmolyte-mediated alterations in interfacial water structures. These observations based on an inorganic interface are reciprocated in amyloid beta 40 dimerization dynamics, highlighting osmolytes' efficacy in mitigating salt-induced aggregation. A molecular analysis suggests that the differential modes of interaction, as found here for glycine and betaine, are prevalent across classes of zwitterionic osmolytes.

Effects of Crumpling Stage and Porosity of Graphene Electrode on the Performance of Electrochemical Supercapacitor

The performance characteristics of supercapacitors composed of crumpled graphene electrodes and aqueous NaCl electrolytes are investigated through Molecular Dynamics (MD) simulations using a newly developed crumpled graphene-based supercapacitor model. Results suggest that the three-dimensional configuration of crumpled graphene boosts electrolyte-electrode interaction. This improved interaction, which includes a larger ion-accessible zone, increases the specific capacitance of the supercapacitor by roughly 400% (16.4 μF/cm2) compared to planar graphene electrodes. Examining the effect of different stages of crumpling and the inclusion of pores on the electrode surface shows that the stages of crumpling substantially influence the supercapacitor performance. A smaller crumpling radius, meaning fully crumpled stage, improves the performance as increased crumpling leads to better packing efficiency, which aids in more ion separation. Furthermore, adding pores on the surface of crumpled graphene improves the ion accessibility by creating additional adsorption sites. An exceptional capacitance of 19.8 μF/cm2 is obtained for a porosity of 20%. However, the results suggest that the in-plane-porosity of the electrode needs to be optimized as there is a decline in specific capacitance after that point (20% porosity), indicating a suboptimal relationship between the charge distribution, specific surface area (SSA) and the porosity of the electrode.

Buffer Screening of Protein Formulations Using a Coarse-Grained Protocol Based on Medicinal Chemistry Interactions

In drug and vaccine development, the designed protein formulation should be highly stable against the temperature, pH, buffer, excipients, and other environmental settings. Similarly, in a sensing unit, one needs to know how strongly two biomolecules bind to guide the design of the biorecognition unit accordingly. Typically, the community performs a series of experiments to thoroughly examine the parameter space, the so-called design-of-experiment (DoE) method, to identify the optimal formulation conditions. Unfortunately, extensive physical testing entails high costs, repeatability issues, and a lack of in-depth knowledge of the underlying mechanisms that affect the final outcome. To address these challenges, we developed a physics-based simulation protocol for buffer screening of protein formulations. We are introducing a coarse-grained molecular simulation protocol that consists of six different interactions. The so-called medicinal chemistry interactions (electrostatics, hydrophobicity, hydrogen bonding propensity, disulfide bonding, and water-water) are based on the physical nature of the protein's amino acid and the partitioning/polarity of any other chemical constituent. The protocol is applied in immunoglobulin-based monoclonal antibodies. We have analyzed the protein behavior as a function of acidity (pH) to discover the isoelectric point by solving the Poisson-Boltzmann equation in a mesoscale grid. To identify the conditions under which the protein oligomerizes in a given buffer, pH, temperature, and ionic strength, we are performing dissipative particle dynamics (DPD) simulations. The protocol allows researchers to reach the high time/space scales required to study protein formulations in their full complexity. Combined with the disruptive protein folding artificial intelligence (AI) algorithms that have been recently developed, the protocol creates a powerful digital framework for cultivating advanced pharmaceutical and biological applications.

Structural Rearrangement of the AT1 Receptor Modulated by Membrane Thickness and Tension

Membrane-embedded mechanosensitive (MS) proteins, including ion channels and G-protein coupled receptors (GPCRs), are essential for the transduction of external mechanical stimuli into biological signals. The angiotensin II type 1 (AT1) receptor plays many important roles in cardiovascular regulation and is associated with diseases such as hypertension and congestive heart failure. The membrane-mediated activation of the AT1 receptor is not well understood, despite this being one of the most widely studied GPCRs within the context of biased agonism. Here, we use extensive molecular dynamics (MD) simulations to characterize the effect of the local membrane environment on the activation of the AT1 receptor. We show that membrane thickness plays an important role in the stability of active and inactive states of the receptor, as well as the dynamic interchange between states. Furthermore, our simulation results show that membrane tension is effective in driving large-scale structural changes in the inactive state such as the outward movement of transmembrane helix 6 to stabilize intermediate active-like conformations. We conclude by comparing our simulation observations with AlphaFold 2 predictions, as a proxy to experimental structures, to provide a framework for how membrane mediated stimuli can facilitate activation of the AT1 receptor through the β-arrestin signaling pathway.

Alzheimer's Disease Immunotherapy and Mimetic Peptide Design for Drug Development: Mutation Screening, Molecular Dynamics, and a Quantum Biochemistry Approach Focusing on Aducanumab::Aβ2-7 Binding Affinity

Seven treatments are approved for Alzheimer's disease, but five of them only relieve symptoms and do not alter the course of the disease. Aducanumab (Adu) and lecanemab are novel disease-modifying antiamyloid-β (Aβ) human monoclonal antibodies that specifically target the pathophysiology of Alzheimer's disease (AD) and were recently approved for its treatment. However, their administration is associated with serious side effects, and their use is limited to early stages of the disease. Therefore, drug discovery remains of great importance in AD research. To gain new insights into the development of novel drugs for Alzheimer's disease, a combination of techniques was employed, including mutation screening, molecular dynamics, and quantum biochemistry. These were used to outline the interfacial interactions of the Aducanumab::Aβ2-7 complex. Our analysis identified critical stabilizing contacts, revealing up to 40% variation in the affinity of the Adu chains for Aβ2-7 depending on the conformation outlined. Remarkably, two complementarity determining regions (CDRs) of the Adu heavy chain (HCDR3 and HCDR2) and one CDR of the Adu light chain (LCDR3) accounted for approximately 77% of the affinity of Adu for Aβ2-7, confirming their critical role in epitope recognition. A single mutation, originally reported to have the potential to increase the affinity of Adu for Aβ2-7, was shown to decrease its structural stability without increasing the overall binding affinity. Mimetic peptides that have the potential to inhibit Aβ aggregation were designed by using computational outcomes. Our results support the use of these peptides as promising drugs with great potential as inhibitors of Aβ aggregation.

Enhancement of tryptophan 2-monooxygenase thermostability by semi-rational enzyme engineering: a strategic design to minimize experimental investigation

Tryptophan 2-monooxygenase (TMO) is an FAD-bound flavoenzyme which catalyzes the oxidative decarboxylation of l-tryptophan to produce indole-3-acetamide (IAM) and carbon dioxide. The reaction of TMO is the first step of indole-3-acetic acid (IAA) biosynthesis. Although TMO is of interest for mechanistic studies and synthetic biology applications, the enzyme has low thermostability and soluble expression yield. Herein, we employed a combined approach of rational design using computational tools with site-saturation mutagenesis to screen for TMO variants with significantly improved thermostability properties and soluble protein expression. The engineered TMO variants, TMO-PWS and TMO-PWSNR, possess melting temperatures (T m) of 65 °C, 17 °C higher than that of the wild-type enzyme (TMO-WT). At 50 °C, the stabilities (t 1/2) of TMO-PWS and TMO-PWSNR were 85-fold and 92.4-fold higher, while their soluble expression yields were 1.4-fold and 2.1-fold greater than TMO-WT, respectively. Remarkably, the kinetic parameters of these variants were similar to those of the wild-type enzymes, illustrating that they are promising candidates for future studies. Molecular dynamic simulations of the wild-type and thermostable TMO variants identified key interactions for enhancing these improvements in the biophysical properties of the TMO variants. The introduced mutations contributed to hydrogen bond formation and an increase in the regional hydrophobicity, thereby, strengthening the TMO structure.

Molecular origins of absorption wavelength variation among phycocyanobilin-binding proteins

Phycocyanobilin (PCB)-binding proteins, including cyanobacteriochromes and phytochromes, function as photoreceptors and exhibit a wide range of absorption maximum wavelengths. To elucidate the color-tuning mechanisms among these proteins, we investigated seven crystal structures of six PCB-binding proteins: Anacy_2551g3, AnPixJg2, phosphorylation-responsive photosensitive histidine kinase, RcaE, Sb.phyB(PG)-PCB, and Slr1393g3. Employing a quantum chemical/molecular mechanical approach combined with a polarizable continuum model, our analysis revealed that differences in absorption wavelengths among PCB-binding proteins primarily arise from variations in the shape of the PCB molecule itself, accounting for a ∼150 nm difference. Remarkably, calculated excitation energies sufficiently reproduced the absorption wavelengths of these proteins spanning ∼200 nm, including 728 nm for Anacy_2551g3. However, assuming the hypothesized lactim conformation resulted in a significant deviation from the experimentally measured absorption wavelength for Anacy_2551g3. The significantly red-shifted absorption wavelength of Anacy_2551g3 can unambiguously be explained by the significant overlap of molecular orbitals between the two pyrrole rings at both edges of the PCB chromophore without the need to hypothesize lactim formation.

Wetting Preference of Silica Surfaces in the Context of Underground Hydrogen Storage: A Molecular Dynamics Perspective

The growing interest in large-scale underground hydrogen (H2) storage (UHS) emphasizes the need for a comprehensive understanding of the fundamental characteristics of subsurface environments. The wetting preference of subsurface rock is a crucial parameter influencing the H2 flow behavior during storage and withdrawal processes. In this study, we utilized molecular dynamics simulation to evaluate the wetting preference of the silica surface in subsurface hydrogen systems, with the aim of addressing disparities observed in experimental results. We conducted an initial comprehensive assessment of potential models, comparing the wettability of five common silica surfaces with different surface morphologies and hydroxyl densities in CO2-H2/water/silica systems against experimental data. After introducing the INTERFACE force field as the most accurate potential model for the silica surface, we evaluated the wetting behavior of the α-quartz (101) surface with a hydroxyl density of 5.9 number/nm2 under the impact of actual geological storage conditions (333-413 K and 10-30 MPa), the coexistence of cushion gases (i.e., CO2, CH4, and N2) at various mole fractions, and pH levels ranging from 2 to 11 characterized through considering the negative charges of 0 to -0.12 C/m2 via deprotonation of silanol on the silica surface. Our results indicate that neither pressure nor temperature has a significant impact on the wetness of the silica in the case of pure H2 (single component UHS operations). However, when CO2 coexists with H2, especially at higher mole fractions, an increase in pressure and a decrease in temperature lead to higher contact angles. Moreover, when the mole fraction of cushion gas ranges from 0 to 1, the contact angle increases 20, 9.5, and 4.5° for CO2, CH4, and N2, respectively, on the neutral silica substrate. Interestingly, at higher pH conditions where the silica surface carries a negative charge, the contact angle considerably reduces where surface charges of -0.03 and -0.06 C/m2 result in an average reduction of 20 and 80% in the contact angle, respectively. More importantly, at a pH of ∼11 (-0.12 C/m2), a 0° contact angle is observed for the silica surface under all temperatures, pressures, types of cushion gases, and varying mole fractions.

Flow-based bioconjugation of coumarin phosphatidylethanolamine probes: Optimised synthesis and membrane molecular dynamics studies

A series of phosphatidylethanolamine fluorescent probes head-labelled with 3-carboxycoumarin was prepared by an improved bioconjugation approach through continuous flow synthesis. The established procedure, supported by a design of experiment (DoE) set-up, resulted in a significant reduction in the reaction time compared to the conventional batch method, in addition to a minor yield increase. The characterization of these probes was enhanced by an in-depth molecular dynamics (MD) study of the behaviour of a representative probe of this family, 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine labelled with 3-carboxycoumarin (POPE-COUM), in bilayers of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC)/1-stearoyl-2-linoleoyl-sn-glycero-3-phosphocholine (SLPC) 2:1, mimicking the composition of the egg yolk lecithin membranes recently used experimentally by our group to study POPE-COUM as a biomarker of the oxidation state and integrity of large unilamellar vesicles (LUVs). The MD simulations revealed that the coumarin group is oriented towards the bilayer interior, leading to a relatively internal location, in agreement with what is observed in the nitrobenzoxadiazole fluorophore of commercial head-labelled NBD-PE probes. This behaviour is consistent with the previously stated hypothesis that POPE-COUM is entirely located within the LUVs structure. Hence, the delay on the oxidation of the probe in the oxygen radical absorbance capacity (ORAC) assays performed is related with the inaccessibility of the probe until alteration of the LUV structure occurs. Furthermore, our simulations show that POPE-COUM exerts very little global and local perturbation on the host bilayer, as evaluated by key properties of the unlabelled lipids. Together, our findings establish PE-COUM as suitable fluorescent lipid analogue probes.

Effect of two activators on the gating of a K2P channel

TWIK-related potassium channel 1 (TREK1), a two-pore-domain mammalian potassium (K+) channel, regulates the resting potential across cell membranes, presenting a promising therapeutic target for neuropathy treatment. The gating of this channel converges in the conformation of the narrowest part of the pore: the selectivity filter (SF). Various hypotheses explain TREK1 gating modulation, including the dynamics of loops connecting the SF with transmembrane helices and the stability of hydrogen bond (HB) networks adjacent to the SF. Recently, two small molecules (Q6F and Q5F) were reported as activators that affect TREK1 by increasing its open probability in single-channel current measurements. Here, using molecular dynamics simulations, we investigate the effect of these ligands on the previously proposed modulation mechanisms of TREK1 gating compared to the apo channel. Our findings reveal that loop dynamics at the upper region of the SF exhibit only a weak correlation with permeation events/nonpermeation periods, whereas the HB network behind the SF appears more correlated. These nonpermeation periods arise from both distinct mechanisms: a C-type inactivation (resulting from dilation at the top of the SF), which has been described previously, and a carbonyl flipping in an SF binding site. We find that, besides the prevention of C-type inactivation in the channel, the ligands increase the probability of permeation by modulating the dynamics of the carbonyl flipping, influenced by a threonine residue at the bottom of the SF. These results offer insights for rational ligand design to optimize the gating modulation of TREK1 and related K+ channels.

A rare case of uncharacterized autoinflammatory disease: Patient carrying variations in NLRP3 and TNFRSF1A genes

Tumor necrosis factor type 1A receptor-associated periodic syndrome (TRAPS) and cryopyrin-associated autoinflammatory syndrome (CAPS) are rare monogenic autoinflammatory diseases (AIDs) mainly caused by pathogenic variations in the TNFRSF1A and NLRP3 genes, respectively. Here, we describe a unique patient presenting with symptoms overlapping both TRAPS and CAPS, without known pathogenic variants in the respective genes. The patient harbored the p.Val200Met variation in NLRP3 and the p.Ser226Cys variation in TNFRSF1A, prompting us to delve deeper into the functional analysis due to conflicting or inconclusive pathogenicity interpretations of the variants across various databases. Molecular dynamics analysis of the p.Val200Met variation in NLRP3 revealed a rigid conformation in the helical domain 2 subdomain of the NACHT domain. This increased rigidity suggests a potential mechanism by which this variation supports the assembly of the NLRP3 inflammasome. Notably, the patient's peripheral mononuclear blood cells demonstrated an elevated IL-1β response upon lipopolysaccharides (LPS) induction. Subsequent initiation of anti-IL-1β therapy resulted in a significant alleviation of the patient's symptoms, further supporting our hypothesis. We interpret these findings as suggestive of a potential pathophysiological role for the NLPR3 p.Val200Met variation in shaping the patient's clinical phenotype, which was also supported by clinical and genetic analysis of the family. This case underscores the complexity of the genetic landscape in AIDs and highlights the value of combining family genetic and functional data to refine the understanding and management of such challenging cases.

Bridging in silico and in vitro perspectives to unravel molecular mechanisms underlying the inhibition of β-glucuronidase by coumarins from Hibiscus trionum

Unraveling the intricacies of β-glucuronidase inhibition is pivotal for developing effective strategies in applications specific to gastrointestinal health and drug metabolism. Our study investigated the efficacy of some Hibiscus trionum phytochemicals as β-glucuronidase inhibitors. The results showed that cleomiscosin A and mansonone H emerged as the most potent inhibitors, with IC50 values of 3.97 ± 0.35 μM and 10.32 ± 1.85 μM, respectively. Mechanistic analysis of β-glucuronidase inhibition indicated that cleomiscosin A and the reference drug EGCG displayed a mixed inhibition mode against β-glucuronidase, while mansonone H exhibited noncompetitive inhibition against β-glucuronidase. Docking studies revealed that cleomiscosin A and mansonone H exhibited the lowest binding affinities, occupying the same site as EGCG, and engaged significant key residues in their binding mechanisms. Using a 30 ns molecular dynamics (MD) simulation, we explored the interaction dynamics of isolated compounds with β-glucuronidase. Analysis of various MD parameters showed that cleomiscosin A and mansonone H exhibited consistent trajectories and significant energy stabilization with β-glucuronidase. These computational insights complemented experimental findings, underscoring the potential of cleomiscosin A and mansonone H as β-glucuronidase inhibitors.

Mica Lattice Orientation of Epitaxially Grown Amyloid β25-35 Fibrils

β-amyloid (Aβ) peptides form self-organizing fibrils in Alzheimer's disease. The biologically active, toxic Aβ25-35 fragment of the full-length Aβ-peptide forms a stable, oriented filament network on the mica surface with an epitaxial mechanism at the timescale of seconds. While many of the structural and dynamic features of the oriented Aβ25-35 fibrils have been investigated before, the β-strand arrangement of the fibrils and their exact orientation with respect to the mica lattice remained unknown. By using high-resolution atomic force microscopy, here, we show that the Aβ25-35 fibrils are oriented along the long diagonal of the oxygen hexagon of mica. To test the structure and stability of the oriented fibrils further, we carried out molecular dynamics simulations on model β-sheets. The models included the mica surface and a single fibril motif built from β-strands. We show that a sheet with parallel β-strands binds to the mica surface with its positively charged groups, but the C-terminals of the strands orient upward. In contrast, the model with antiparallel strands preserves its parallel orientation with the surface in the molecular dynamics simulation, suggesting that this model describes the first β-sheet layer of the mica-bound Aβ25-35 fibrils well. These results pave the way toward nanotechnological construction and applications for the designed amyloid peptides.

Binding of Inhibitors to Nuclear Localization Signal Peptide from Venezuelan Equine Encephalitis Virus Capsid Protein Explored with All-Atom Replica Exchange Molecular Dynamics

Several small molecule inhibitors have been designed to block binding of the Venezuelan equine encephalitis virus (VEEV) nuclear localization signal (NLS) sequence to the importin-α nuclear transport protein. To probe the inhibition mechanism on a molecular level, we used all-atom explicit water replica exchange molecular dynamics to study the binding of two inhibitors, I1 and I2, to the coreNLS peptide, representing the core fragment of the VEEV NLS sequence. Our objective was to evaluate the possibility of masking wherein binding of these inhibitors to the coreNLS occurs prior to its binding to importin-α. We found that the free energy of I1 and I2 binding to the coreNLS is less favorable than that to importin-α. This outcome argues against preemptive inhibitor binding to the coreNLS prior to importin-α. Instead, both inhibitors are expected to compete with the coreNLS peptide for binding to importin-α. The two factors responsible for the low affinities of the inhibitors to the coreNLS peptide are (i) the low cooperativity of binding to the peptide and (ii) the strong hydrophobic effect associated with binding to importin-α. Our results further show that upon binding to the coreNLS peptide, the inhibitors form multiple diverse binding poses. The coreNLS peptide coincubated with I1 and I2 adopts several conformational states, including open and collapsed, which underscores the fluidity of the coreNLS conformational ensemble as a target for inhibitors. Taken together with our prior investigations, this study sheds light on the molecular mechanism by which I1 and I2 ligands inhibit binding of the VEEV capsid protein to importin-α.

Computational and experimental identification of an exceptionally efficient ethyl ester synthetase with broad substrate specificity and high product yield, suggests potential for industrial biocatalysis

Transesterification plays a crucial role in the synthesis of diverse esters in organic synthesis but is barely reported in biocatalysis. Here, we computationally identify the salicylic acid-binding protease 2 (SABP2) as an efficient ethyl ester bond synthetase by QM/MM MD and free energy simulations and present the practical and effective utilization of SABP2 as an eco-friendly biocatalyst for transesterification reactions by a series of experiments. Our findings demonstrate that SABP2 efficiently catalyzes the transesterification reaction between the carboxyl acid group of promiscuous aromatic substrates and ethanol to produce the corresponding ethyl esters. Notably, while SABP2 exhibits its native capability to catalyze the hydrolysis of the methyl salicylate (MeSA), the transesterification rate (producing ethyl salicylate, EtSA) is about 3500 times higher than the hydrolysis rate. Additionally, a range of aromatic methyl esters are employed in the transesterification process, resulting in high yields (up to 98.9 %) of the corresponding ethyl esters. These results indicate a broad substrate scope for SABP2-catalyzed transesterification reactions, demonstrating its potential as a valuable biocatalyst for ester synthesis in industry.

Ligand-Dependent and G Protein-Dependent Properties for the Sweet Taste Heterodimer, TAS1R2/1R3

The heterodimeric sweet taste receptor, TAS1R2/1R3, is a class C G protein-coupled receptor (GPCR) that couples to gustducin (Gt), a G protein (GP) specifically involved in taste processing. This makes TAS1R2/1R3 a possible target for newly developing low caloric ligands that taste sweet to address obesity and diabetes. The activation of TAS1R2/1R3 involves the insertion of the GαP C-terminus of the GP into the GPCR in response to ligand binding. However, it is not known for sure whether the GP inserts into the TAS1R2 or TAS1R3 intracellular region of this GPCR dimer. Moreover, TAS1R2/1R3 can also connect to other GPs, such as Gs, Gi1, Gt3, Go, Gq, and G12. These GPs have different C-termini that may modify GPCR signaling. To understand the possible GP dependence of sweet perception, we use molecular dynamic (MD) simulations to examine the coupling of various GαP C20 termini to TAS1R2/1R3 for various steviol glycoside ligands and an artificial sweetener. Since the C20 could interact with the transmembrane domain (TMD) of either TAS1R2 (TMD2) or TAS1R3 (TMD3), we consider both cases. Without any sweetener, we find that the apo GPCR shows similar Go and Gt selectivities, while all steviol glycoside ligands increase the selectivity of Gt but decrease Go selectivity at TMD2. Interestingly, we find that high sweet rebaudioside M (RebM) and RebD ligands show better interactions of C20 at TMD3 for the Gt protein, but low sweet RebC and hydRebM ligands show better interaction of C20 at TMD2 for the Gt protein. Thus, our MD simulation suggests that TAS1R2/1R3 may couple the GP to either 1R2 or to 1R3 and that it can couple other GPs compared to Gt. This will likely lead to multimodal functions producing multiple patterns of intracellular signaling for sweet taste receptors, depending on the particular sweetener. Directing the GP to one of the other may have beneficial therapeutic outcomes.

Lipophilic compounds restore function to neurodevelopmental-associated KCNQ3 mutations

A major driver of neuronal hyperexcitability is dysfunction of K+ channels, including voltage-gated KCNQ2/3 channels. Their hyperpolarized midpoint of activation and slow activation and deactivation kinetics produce a current that regulates membrane potential and impedes repetitive firing. Inherited mutations in KCNQ2 and KCNQ3 are linked to a wide spectrum of neurodevelopmental disorders (NDDs), ranging from benign familial neonatal seizures to severe epileptic encephalopathies and autism spectrum disorders. However, the impact of these variants on the molecular mechanisms underlying KCNQ3 channel function remains poorly understood and existing treatments have significant side effects. Here, we use voltage clamp fluorometry, molecular dynamic simulations, and electrophysiology to investigate NDD-associated variants in KCNQ3 channels. We identified two distinctive mechanisms by which loss- and gain-of function NDD-associated mutations in KCNQ3 affect channel gating: one directly affects S4 movement while the other changes S4-to-pore coupling. MD simulations and electrophysiology revealed that polyunsaturated fatty acids (PUFAs) primarily target the voltage-sensing domain in its activated conformation and form a weaker interaction with the channel's pore. Consistently, two such compounds yielded partial and complete functional restoration in R227Q- and R236C-containing channels, respectively. Our results reveal the potential of PUFAs to be developed into therapies for diverse KCNQ3-based channelopathies.

The retinal chromophore environment in an inward light-driven proton pump studied by solid-state NMR and hydrogen-bond network analysis

Inward proton pumping is a relatively new function for microbial rhodopsins, retinal-binding light-driven membrane proteins. So far, it has been demonstrated for two unrelated subgroups of microbial rhodopsins, xenorhodopsins and schizorhodopsins. A number of recent studies suggest unique retinal-protein interactions as being responsible for the reversed direction of proton transport in the latter group. Here, we use solid-state NMR to analyze the retinal chromophore environment and configuration in an inward proton-pumping Antarctic schizorhodopsin. Using fully 13C-labeled retinal, we have assigned chemical shifts for every carbon atom and, assisted by structure modelling and molecular dynamics simulations, made a comparison with well-studied outward proton pumps, identifying locations of the unique protein-chromophore interactions for this functional subclass of microbial rhodopsins. Both the NMR results and molecular dynamics simulations point to the distinctive polar environment in the proximal part of the retinal, which may result in a hydration pattern dramatically different from that of the outward proton pumps, causing the reversed proton transport.

Computational Insights into the Non-Heme Diiron Alkane Monooxygenase Enzyme AlkB: Electronic Structures, Dioxygen Activation, and Hydroxylation Mechanism of Liquid Alkanes

Alkane monooxygenase (AlkB) is a membrane-spanning metalloenzyme that catalyzes the terminal hydroxylation of straight-chain alkanes involved in the microbially mediated degradation of liquid alkanes. According to the cryoEM structures, AlkB features a unique multihistidine ligand coordination environment with a long Fe-Fe distance in its active center. Up to now, how AlkB employs the diiron center to activate dioxygen and which species is responsible for triggering the hydroxylation are still elusive. In this work, we constructed computational models and performed quantum mechanics/molecular mechanics (QM/MM) calculations to illuminate the electronic characteristics of the diiron active center and how AlkB carries out the terminal hydroxylation. Our calculations revealed that the spin-spin interaction between two irons is rather weak. The dioxygen may ligate to either the Fe1 or Fe2 atom and prefers to act as a linker to increase the spin-spin interaction of two irons, facilitating the dioxygen cleavage to generate the highly reactive Fe(IV)═O. Thus, AlkB employs Fe(IV)═O to trigger the hydrogen abstraction. In addition, the previously suggested mechanism that AlkB uses both the dioxygen and Fe-coordinated water to perform hydroxylation was calculated to be unlikely. Besides, our results indicate that AlkB cannot use the Fe-coordinated dioxygen to directly trigger hydrogen abstraction.

Cholesterol inhibits oxygen permeation through biological membranes: mechanism against double-bond peroxidation

The presence of oxygen molecules (O2) in biological membranes promotes lipid peroxidation of phospholipids with unsaturated acyl chains. On the other hand, cholesterol is considered to be an antioxidant molecule as it has a significant barrier effect on the permeation of O2 across membranes. However, a comprehensive explanation of how cholesterol affects the distribution and diffusion of O2 within lipid bilayers is yet to be established. In this study, we investigated the interaction of oxygen molecules with polyunsaturated lipid bilayers using molecular dynamics (MD) simulations. The degree of lipid unsaturation and the concentration of cholesterol were varied to study the permeation of O2. The free energy profile of O2 diffusing from the water phase to the lipid bilayer was calculated using biased umbrella MD simulations. The results show that O2 passively translocates into the membrane without changing the physical properties of the bilayer. Interestingly, in the unsaturated lipid bilayers the presence of cholesterol led to a significantly decreased permeation of O2 and an increase in the lipid chain order. Our results indicate that the hydroxyl groups of cholesterol strongly interact with the O2 molecules effectively inhibiting interactions between the oxygens and the double bonds in unsaturated lipid tails. In addition, a linear relationship between permeation and the ratio of membrane thickness and area per lipid was found. These insights can help our understanding of how the degree of unsaturation in a lipid tail and cholesterol affect lipid peroxidation at the molecular level.

Exploring Protonation State, Ion Binding, and Photoactivated Channel Opening of an Anion Channelrhodopsin by Molecular Simulations

Channelrhodopsins are light-gated ion channels with a retinal chromophore found in microbes and are widely used in optogenetics, a field of neuroscience that utilizes light to regulate neuronal activity. GtACR1, an anion conducting channelrhodopsin derived from Guillardia theta, has attracted attention for its application as a neuronal silencer in optogenetics because of its high conductivity and selectivity. However, atomistic mechanisms of channel photoactivation and ion conduction have not yet been elucidated. In the present study, we investigated the molecular characteristics of GtACR1 and its photoactivation processes by molecular simulations. The QM/MM RWFE-SCF method which combines highly accurate quantum chemistry calculations with long-time molecular dynamics (MD) simulations were used to model protein structures of the wild-type and mutants with different protonation states of key groups and to calculate absorption energies for verification of the models. The QM/MM modeling together with MD simulations of free-energy calculations favors protonation of a key counterion carboxyl group of Asp234 with a strong binding of a chloride ion in the extracellular pocket in the dark state. A channel open state was also successfully modeled by the QM/MM RWFE-SCF free-energy optimizations, providing atomistic insights into the channel activation mechanism.

Statistical accuracy of molecular dynamics-based methods for sampling conformational ensembles of disordered proteins

The characterization of the statistical ensemble of conformations of intrinsically disordered regions (IDRs) is a great challenge both from experimental and computational points of view. In this respect, a number of protocols have been developed using molecular dynamics (MD) simulations to sample the huge conformational space of the molecule. In this work, we consider one of the best methods available, replica exchange solute tempering (REST), as a reference to compare the results obtained using this method with the results obtained using other methods, in terms of experimentally measurable quantities. Along with the methods assessed, we propose here a novel protocol called probabilistic MD chain growth (PMD-CG), which combines the flexible-meccano and hierarchical chain growth methods with the statistical data obtained from tripeptide MD trajectories as the starting point. The system chosen for testing is a 20-residue region from the C-terminal domain of the p53 tumor suppressor protein (p53-CTD). Our results show that PMD-CG provides an ensemble of conformations extremely quickly, after suitable computation of the conformational pool for all peptide triplets of the IDR sequence. The measurable quantities computed on the ensemble of conformations agree well with those based on the REST conformational ensemble.

Development of an anisotropic polarizable model for the all-atom AMOEBA force field

For planar and rigid π-conjugated molecular systems, electrostatic and inductive interactions are pivotal in governing molecular packing structures and electron polarization energies. These electrostatic interactions typically exhibit an anisotropic nature within π-conjugated systems. In this study, we utilize the atoms in molecules (AIM) theory in conjunction with linear response theory to decompose molecular polarizability into distributed atomic polarizability tensors. On the basis of atomic polarizability tensors, we extended an anisotropic polarizable model into the AMOEBA polarizable force field. Both anisotropic and isotropic polarizable models in combination with various density functional theory (DFT)-derived atomic multipoles were applied to optimize the experimental crystals of naphthalene and anthracene. Furthermore, these two types of electrostatic models, coupled with the evolutionary algorithm USPEX program, are utilized to predict the crystal structures of oligoacenes. Our findings demonstrate that the anisotropic polarizable model exhibits superior performance in crystal refinement and crystal structure prediction. This enriched anisotropic polarizable model is seamlessly integrated into the AMOEBA polarizable force field and readily applicable within our modified Tinker program.

van der Waals Radii of Free and Bonded Atoms from Hydrogen (Z = 1) to Oganesson (Z = 118)

Reliable numerical values of van der Waals (vdW) radii are required for constructing empirical force fields, vdW-inclusive density functional, and quantum-chemical methods, as well as for implicit solvent models. However, multiple definitions exist for vdW radii, involving either equilibrium or the closest contact distances between free or bonded atoms within molecules or crystals. For the paradigmatic case of the hydrogen atom, its reported vdW radius fluctuates between 2.15 and 3.70 Bohr depending on the definition, leading to a high uncertainty in calculations and different conceptual interpretations of noncovalent interactions. In this work, we systematically review different definitions and methodologies to establish the free and bonded vdW radii for hydrogen, based on equilibrium vdW distances in noncovalently bonded molecules, enveloping electron density cutoffs, noncovalent positron bonds in hydrogen anion dimer, vacuum virtual photon cloud caused by the hydrogen atom, and atomic dipole polarizability. By doing so, we show that the vdW radius of the free hydrogen atom is 3.16 ± 0.06 Bohr. By employing the most general and elegant definition of atomic vdW radius as a function of the atomic polarizability, we tabulate consistent values of vdW radii for all atoms in the periodic table up to Z = 118.

Effective Inclusion of Electronic Polarization Improves the Description of Electrostatic Interactions: The prosECCo75 Biomolecular Force Field

prosECCo75 is an optimized force field effectively incorporating electronic polarization via charge scaling. It aims to enhance the accuracy of nominally nonpolarizable molecular dynamics simulations for interactions in biologically relevant systems involving water, ions, proteins, lipids, and saccharides. Recognizing the inherent limitations of nonpolarizable force fields in precisely modeling electrostatic interactions essential for various biological processes, we mitigate these shortcomings by accounting for electronic polarizability in a physically rigorous mean-field way that does not add to computational costs. With this scaling of (both integer and partial) charges within the CHARMM36 framework, prosECCo75 addresses overbinding artifacts. This improves agreement with experimental ion binding data across a broad spectrum of systems─lipid membranes, proteins (including peptides and amino acids), and saccharides─without compromising their biomolecular structures. prosECCo75 thus emerges as a computationally efficient tool providing enhanced accuracy and broader applicability in simulating the complex interplay of interactions between ions and biomolecules, pivotal for improving our understanding of many biological processes.

Molecular Basis of the Substrate Specificity of Phosphotriesterase from Pseudomonas diminuta: A Combined QM/MM MD and Electron Density Study

The occurrence of organophosphorus compounds, pesticides, and flame-retardants in wastes is an emerging ecological problem. Bacterial phosphotriesterases are capable of hydrolyzing some of them. We utilize modern molecular modeling tools to study the hydrolysis mechanism of organophosphorus compounds with good and poor leaving groups by phosphotriesterase from Pseudomonas diminuta (Pd-PTE). We compute Gibbs energy profiles for enzymes with different cations in the active site: native Zn2+cations and Co2+cations, which increase the steady-state rate constant. Hydrolysis occurs in two elementary steps via an associative mechanism and formation of the pentacoordinated intermediate. The first step, a nucleophilic attack, occurs with a low energy barrier independently of the substrate. The second step has a low energy barrier and considerable stabilization of products for substrates with good leaving groups. For substrates with poor leaving groups, the reaction products are destabilized relative to the ES complex that suppresses the reaction. The reaction proceeds with low energy barriers for substrates with good leaving groups with both Zn2+and Co2+cations in the active site; thus, the product release is likely to be a limiting step. Electron density and geometry analysis of the QM/MM MD trajectories of the intermediate states with all considered compounds allow us to discriminate substrates by their ability to be hydrolyzed by the Pd-PTE. For hydrolyzable substrates, the cleaving bond between a phosphorus atom and a leaving group is elongated, and electron density depletion is observed on the Laplacian of electron density maps.