OpenMM 4: A Reusable, Extensible, Hardware Independent Library for High Performance Molecular Simulation

Plugin-based interoperability and ecosystem management for the MolSSI Driver Interface Project

The MolSSI Driver Interface (MDI) Project is an effort to simplify and standardize the process of enabling tight interoperability between independently developed code bases and is supported by numerous software packages across the domain of chemical physics. It enables a wide variety of use cases, including quantum mechanics/molecular mechanics, advanced sampling, path integral molecular dynamics, machine learning, ab initio molecular dynamics, etc. We describe two major developments within the MDI Project that provide novel solutions to key interoperability challenges. The first of these is the development of the MDI Plugin System, which allows MDI-supporting libraries to be used as highly modular plugins, with MDI enforcing a standardized application programming interface across plugins. Codes can use these plugins without linking against them during their build process, and end-users can select which plugin(s) they wish to use at runtime. The MDI Plugin System features a sophisticated callback system that allows codes to interact with plugins on a highly granular level and represents a significant advancement toward increased modularity among scientific codes. The second major development is MDI Mechanic, an ecosystem management tool that utilizes Docker containerization to simplify the process of developing, validating, maintaining, and deploying MDI-supporting codes. Additionally, MDI Mechanic provides a framework for launching MDI simulations in which each interoperating code is executed within a separate computational environment. This eliminates the need to compile multiple production codes within a single computational environment, reducing opportunities for dependency conflicts and lowering the barrier to entry for users of MDI-enabled codes.

PyRETIS 3: Conquering rare and slow events without boundaries

We present and discuss the advancements made in PyRETIS 3, the third instalment of our Python library for an efficient and user-friendly rare event simulation, focused to execute molecular simulations with replica exchange transition interface sampling (RETIS) and its variations. Apart from a general rewiring of the internal code towards a more modular structure, several recently developed sampling strategies have been implemented. These include recently developed Monte Carlo moves to increase path decorrelation and convergence rate, and new ensemble definitions to handle the challenges of long-lived metastable states and transitions with unbounded reactant and product states. Additionally, the post-analysis software PyVisa is now embedded in the main code, allowing fast use of machine-learning algorithms for clustering and visualising collective variables in the simulation data.

Distinct roles of the major binding residues in the cation-binding pocket of the melibiose transporter MelB

Salmonella enterica serovar Typhimurium melibiose permease (MelBSt) is a prototype of the major facilitator superfamily (MFS) transporters, which play important roles in human health and diseases. MelBSt catalyzed the symport of galactosides with either Na+, Li+, or H+, but prefers the coupling with Na+. Previously, we determined the structures of the inward- and outward-facing conformation of MelBSt and the molecular recognition for galactoside and Na+. However, the molecular mechanisms for H+- and Na+-coupled symport still remain poorly understood. In this study, we solved two x-ray crystal structures of MelBSt, the cation-binding site mutants D59C at an unliganded apo-state and D55C at a ligand-bound state, and both structures display the outward-facing conformations virtually identical as published. We determined the energetic contributions of three major Na+-binding residues for the selection of Na+ and H+ by free energy simulations. Transport assays showed that the D55C mutant converted MelBSt to a solely H+-coupled symporter, and together with the free-energy perturbation calculation, Asp59 is affirmed to be the sole protonation site of MelBSt. Unexpectedly, the H+-coupled melibiose transport exhibited poor activities at greater bulky ΔpH and better activities at reversal ΔpH, supporting the novel theory of transmembrane-electrostatically localized protons and the associated membrane potential as the primary driving force for the H+-coupled symport mediated by MelBSt. This integrated study of crystal structure, bioenergetics, and free energy simulations, demonstrated the distinct roles of the major binding residues in the cation-binding pocket of MelBSt.

Structure prediction of protein-ligand complexes from sequence information with Umol

Protein-ligand docking is an established tool in drug discovery and development to narrow down potential therapeutics for experimental testing. However, a high-quality protein structure is required and often the protein is treated as fully or partially rigid. Here we develop an AI system that can predict the fully flexible all-atom structure of protein-ligand complexes directly from sequence information. We find that classical docking methods are still superior, but depend upon having crystal structures of the target protein. In addition to predicting flexible all-atom structures, predicted confidence metrics (plDDT) can be used to select accurate predictions as well as to distinguish between strong and weak binders. The advances presented here suggest that the goal of AI-based drug discovery is one step closer, but there is still a way to go to grasp the complexity of protein-ligand interactions fully. Umol is available at: https://github.com/patrickbryant1/Umol .

NPEX: Never give up protein exploration with deep reinforcement learning

Elucidating unknown structures of proteins, such as metastable states, is critical in designing therapeutic agents. Protein structure exploration has been performed using advanced computational methods, especially molecular dynamics and Markov chain Monte Carlo simulations, which require untenably long calculation times and prior structural knowledge. Here, we developed an innovative method for protein structure determination called never give up protein exploration (NPEX) with deep reinforcement learning. The NPEX method leverages the soft actor-critic algorithm and the intrinsic reward system, effectively adding a bias potential without the need for prior knowledge. To demonstrate the method's effectiveness, we applied it to four models: a double well, a triple well, the alanine dipeptide, and the tryptophan cage. Compared with Markov chain Monte Carlo simulations, NPEX had markedly greater sampling efficiency. The significantly enhanced computational efficiency and lack of prior domain knowledge requirements of the NPEX method will revolutionize protein structure exploration.

Design and synthesis of new dihydropyrimidine/sulphonamide hybrids as promising anti-inflammatory agents via dual mPGES-1/5-LOX inhibition

A novel series of dihydropyrimidine/sulphonamide hybrids 3a-j with anti-inflammatory properties have been developed and tested as dual mPGES-1/5-LOX inhibitors. In vitro assay, results showed that compounds 3c, 3e, 3h, and 3j were the most effective dual inhibitors of mPGES-1 and 5-LOX activities. Compound 3j was the most potent dual inhibitor with IC50 values of 0.92 µM and 1.98 µM, respectively. In vivo, anti-inflammatory studies demonstrated that compounds 3c, 3e, 3h, and 3e had considerable anti-inflammatory activity, with EI% ranging from 29% to 71%. Compounds 3e and 3j were equivalent to celecoxib after the first hour but exhibited stronger anti-inflammatory effects than celecoxib after the third and fifth hours. Moreover, compounds 3e and 3j significantly reduced the levels of pro-inflammatory cytokines (PGE2, TNF-α, and IL-6) with gastrointestinal safety profiles. Molecular docking simulations explored the most potent derivatives' binding affinities and interaction patterns within mPGES-1 and 5-LOX active sites. This study disclosed that compound 3j is a promising anti-inflammatory lead with dual mPGES-1/5-LOX inhibition that deserves further preclinical investigation.

Generalized Helmholtz model describes capacitance profiles of ionic liquids and concentrated aqueous electrolytes

In this work, we propose and validate a generalization of the Helmholtz model that can account for both "bell-shaped" and "camel-shaped" differential capacitance profiles of concentrated electrolytes, the latter being characteristic of ionic liquids. The generalization is based on introducing voltage dependence of both the dielectric constant "ϵr(V)" and thickness "L(V)" of the inner Helmholtz layer, as validated by molecular dynamics (MD) simulations. We utilize MD simulations to study the capacitance profiles of three different electrochemical interfaces: (1) graphite/[BMIm+][BF4-] ionic liquid interface; (2) Au(100)/[BMIm+][BF4-] ionic liquid interface; (3) Au(100)/1M [Na+][Cl-] aqueous interface. We compute the voltage dependence of ϵr(V) and L(V) and demonstrate that the generalized Helmholtz model qualitatively describes both camel-shaped and bell-shaped differential capacitance profiles of ionic liquids and concentrated aqueous electrolytes (in lieu of specific ion adsorption). In particular, the camel-shaped capacitance profile that is characteristic of ionic liquid electrolytes arises simply from combination of the voltage-dependent trends of ϵr(V) and L(V). Furthermore, explicit analysis of the inner layer charge density for both concentrated aqueous and ionic liquid double layers reveal similarities, with these charge distributions typically exhibiting a dipolar region closest to the electrode followed by a monopolar peak at larger distances. It is appealing that a generalized Helmholtz model can provide a unified description of the inner layer structure and capacitance profile for seemingly disparate aqueous and ionic liquid electrolytes.

Potential energy landscape of a flexible water model: Equation of state, configurational entropy, and Adam-Gibbs relationship

The potential energy landscape (PEL) formalism is a tool within statistical mechanics that has been used in the past to calculate the equation of states (EOS) of classical rigid model liquids at low temperatures, where computer simulations may be challenging. In this work, we use classical molecular dynamics (MD) simulations and the PEL formalism to calculate the EOS of the flexible q-TIP4P/F water model. This model exhibits a liquid-liquid critical point (LLCP) in the supercooled regime, at (Pc = 150 MPa, Tc = 190 K, and ρc = 1.04 g/cm3) [using the reaction field technique]. The PEL-EOS of q-TIP4P/F water and the corresponding location of the LLCP are in very good agreement with the MD simulations. We show that the PEL of q-TIP4P/F water is Gaussian, which allows us to calculate the configurational entropy of the system, Sconf. The Sconf of q-TIP4P/F water is surprisingly similar to that reported previously for rigid water models, suggesting that intramolecular flexibility does not necessarily add roughness to the PEL. We also show that the Adam-Gibbs relation, which relates the diffusion coefficient D with Sconf, holds for the flexible q-TIP4P/F water model. Overall, our results indicate that the PEL formalism can be used to study molecular systems that include molecular flexibility, the common case in standard force fields. This is not trivial since the introduction of large bending/stretching mode frequencies is problematic in classical statistical mechanics. For example, as shown previously, we find that such high frequencies lead to unphysical (negative) entropy for q-TIP4P/F water when using classical statistical mechanics (yet, the PEL formalism can be applied successfully).

Massively scalable workflows for quantum chemistry: BigChem and ChemCloud

Electronic structure theory, i.e., quantum chemistry, is the fundamental building block for many problems in computational chemistry. We present a new distributed computing framework (BigChem), which allows for an efficient solution of many quantum chemistry problems in parallel. BigChem is designed to be easily composable and leverages industry-standard middleware (e.g., Celery, RabbitMQ, and Redis) for distributed approaches to large scale problems. BigChem can harness any collection of worker nodes, including ones on cloud providers (such as AWS or Azure), local clusters, or supercomputer centers (and any mixture of these). BigChem builds upon MolSSI packages, such as QCEngine to standardize the operation of numerous computational chemistry programs, demonstrated here with Psi4, xtb, geomeTRIC, and TeraChem. BigChem delivers full utilization of compute resources at scale, offers a programable canvas for designing sophisticated quantum chemistry workflows, and is fault tolerant to node failures and network disruptions. We demonstrate linear scalability of BigChem running computational chemistry workloads on up to 125 GPUs. Finally, we present ChemCloud, a web API to BigChem and successor to TeraChem Cloud. ChemCloud delivers scalable and secure access to BigChem over the Internet.

New 2-oxoindole derivatives as multiple PDGFRα/ß and VEGFR-2 tyrosine kinase inhibitors

Two new series of N-aryl acetamides 6a-o and benzyloxy benzylidenes 9a-p based 2-oxoindole derivatives were designed as potent antiproliferative multiple kinase inhibitors. The results of one-dose NCI antiproliferative screening for compounds 6a-o and 9a-p elucidated that the most promising antiproliferative scaffolds were 6f and 9f, which underwent five-dose testing. Notably, the amido congener 6f was the most potent derivative towards pancreatic ductal adenocarcinoma MDA-PATC53 and PL45 cell lines (IC50 = 1.73 µM and 2.40 µM, respectively), and the benzyloxy derivative 9f was the next potent one with IC50 values of 2.85 µM and 2.96 µM, respectively. Both compounds 6f and 9f demonstrated a favorable safety profile when tested against normal prostate epithelial cells (RWPE-1). Additionally, compound 6f displayed exceptional selectivity as a multiple kinase inhibitor, particularly targeting PDGFRα, PDGFRβ, and VEGFR-2 kinases, with IC50 values of 7.41 nM, 6.18 nM, and 7.49 nM, respectively. In contrast, the reference compound Sunitinib exhibited IC50 values of 43.88 nM, 2.13 nM, and 78.46 nM against the same kinases. The derivative 9f followed closely, with IC50 values of 9.9 nM, 6.62 nM, and 22.21 nM for the respective kinases. Both 6f and 9f disrupt the G2/M cell cycle transition by upregulating p21 and reducing CDK1 and cyclin B1 mRNA levels. The interplay between targeted kinases and these cell cycle regulators underpins the G2/M cell cycle arrest induced by our compounds. Also, compounds 6f and 9f fundamentally resulted in entering MDA-PATC53 cells into the early stage of apoptosis with good percentages compared to the positive control Sunitinib. The in silico molecular-docking outcomes of scaffolds 6a-o and 9a-p in VEGFR-2, PDGFRα, and PDGFRβ active sites depicted their ability to adopt essential binding interactions like the reference Sunitinib. Our designed analogs, specifically 6f and 9f, possess promising antiproliferative and kinase inhibitory properties, making them potential candidates for further therapeutic development.

Benchmarking ANI potentials as a rescoring function and screening FDA drugs for SARS-CoV-2 Mpro

Here, we introduce the use of ANI-ML potentials as a rescoring function in the host-guest interaction in molecular docking. Our results show that the "docking power" of ANI potentials can compete with the current scoring functions at the same level of computational cost. Benchmarking studies on CASF-2016 dataset showed that ANI is ranked in the top 5 scoring functions among the other 34 tested. In particular, the ANI predicted interaction energies when used in conjunction with GOLD-PLP scoring function can boost the top ranked solution to be the closest to the x-ray structure. Rapid and accurate calculation of interaction energies between ligand and protein also enables screening of millions of drug candidates/docking poses. Using a unique protocol in which docking by GOLD-PLP, rescoring by ANI-ML potentials and extensive MD simulations along with end state free energy methods are combined, we have screened FDA approved drugs against the SARS-CoV-2 main protease (Mpro). The top six drug molecules suggested by the consensus of these free energy methods have already been in clinical trials or proposed as potential drug molecules in previous theoretical and experimental studies, approving the validity and the power of accuracy in our screening method.

Distinct roles of the major binding residues in the cation-binding pocket of MelB

Salmonella enterica serovar Typhimurium melibiose permease (MelBSt) is a prototype of the major facilitator superfamily (MFS) transporters, which play important roles in human health and diseases. MelBSt catalyzed the symport of galactosides with either H+, Li+, or Na+, but prefers the coupling with Na+. Previously, we determined the structures of the inward- and outward-facing conformation of MelBSt, as well as the molecular recognition for galactoside and Na+. However, the molecular mechanisms for H+- and Na+-coupled symport still remain poorly understood. We have solved two x-ray crystal structures of MelBSt cation-binding site mutants D59C at an unliganded apo-state and D55C at a ligand-bound state, and both structures display the outward-facing conformations virtually identical as published previously. We determined the energetic contributions of three major Na+-binding residues in cation selectivity for Na+ and H+ by the free energy simulations. The D55C mutant converted MelBSt to a solely H+-coupled symporter, and together with the free-energy perturbation calculation, Asp59 is affirmed to be the sole protonation site of MelBSt. Unexpectedly, the H+-coupled melibiose transport with poor activities at higher ΔpH and better activities at reversal ΔpH was observed, supporting that the membrane potential is the primary driving force for the H+-coupled symport mediated by MelBSt. This integrated study of crystal structure, bioenergetics, and free energy simulations, demonstrated the distinct roles of the major binding residues in the cation-binding pocket.

Solid-state NMR MAS CryoProbe enables structural studies of human blood protein vitronectin bound to hydroxyapatite

The low sensitivity of nuclear magnetic resonance (NMR) is a major bottleneck for studying biomolecular structures of complex biomolecular assemblies. Cryogenically cooled probe technology overcomes the sensitivity limitations enabling NMR applications to challenging biomolecular systems. Here we describe solid-state NMR studies of the human blood protein vitronectin (Vn) bound to hydroxyapatite (HAP), the mineralized form of calcium phosphate, using a CryoProbe designed for magic angle spinning (MAS) experiments. Vn is a major blood protein that regulates many different physiological and pathological processes. The high sensitivity of the CryoProbe enabled us to acquire three-dimensional solid-state NMR spectra for sequential assignment and characterization of site-specific water-protein interactions that provide initial insights into the organization of the Vn-HAP complex. Vn associates with HAP in various pathological settings, including macular degeneration eyes and Alzheimer's disease brains. The ability to probe these assemblies at atomic detail paves the way for understanding their formation.

Backbone interactions and secondary structures in phase separation of disordered proteins

Intrinsically disordered proteins (IDPs) are one of the major drivers behind the formation and characteristics of biomolecular condensates. Due to their inherent flexibility, the backbones of IDPs are significantly exposed, rendering them highly influential and susceptible to biomolecular phase separation. In densely packed condensates, exposed backbones have a heightened capacity to interact with neighboring protein chains, which might lead to strong coupling between the secondary structures and phase separation and further modulate the subsequent transitions of the condensates, such as aging and fibrillization. In this mini-review, we provide an overview of backbone-mediated interactions and secondary structures within biomolecular condensates to underscore the importance of protein backbones in phase separation. We further focus on recent advances in experimental techniques and molecular dynamics simulation methods for probing and exploring the roles of backbone interactions and secondary structures in biomolecular phase separation involving IDPs.

Mobile barrier mechanisms for Na+-coupled symport in an MFS sugar transporter

While many 3D structures of cation-coupled transporters have been determined, the mechanistic details governing the obligatory coupling and functional regulations still remain elusive. The bacterial melibiose transporter (MelB) is a prototype of major facilitator superfamily transporters. With a conformation-selective nanobody, we determined a low-sugar affinity inward-facing Na+-bound cryoEM structure. The available outward-facing sugar-bound structures showed that the N- and C-terminal residues of the inner barrier contribute to the sugar selectivity. The inward-open conformation shows that the sugar selectivity pocket is also broken when the inner barrier is broken. Isothermal titration calorimetry measurements revealed that this inward-facing conformation trapped by this nanobody exhibited a greatly decreased sugar-binding affinity, suggesting the mechanisms for substrate intracellular release and accumulation. While the inner/outer barrier shift directly regulates the sugar-binding affinity, it has little or no effect on the cation binding, which is supported by molecular dynamics simulations. Furthermore, the hydron/deuterium exchange mass spectrometry analyses allowed us to identify dynamic regions; some regions are involved in the functionally important inner barrier-specific salt-bridge network, which indicates their critical roles in the barrier switching mechanisms for transport. These complementary results provided structural and dynamic insights into the mobile barrier mechanism for cation-coupled symport.

A continuum of amorphous ices between low-density and high-density amorphous ice

Amorphous ices are usually classified as belonging to low-density or high-density amorphous ice (LDA and HDA) with densities ρLDA ≈ 0.94 g/cm3 and ρHDA ≈ 1.15-1.17 g/cm3. However, a recent experiment crushing hexagonal ice (ball-milling) produced a medium-density amorphous ice (MDA, ρMDA ≈ 1.06 g/cm3) adding complexity to our understanding of amorphous ice and the phase diagram of supercooled water. Motivated by the discovery of MDA, we perform computer simulations where amorphous ices are produced by isobaric cooling and isothermal compression/decompression. Our results show that, depending on the pressure employed, isobaric cooling can generate a continuum of amorphous ices with densities that expand in between those of LDA and HDA (briefly, intermediate amorphous ices, IA). In particular, the IA generated at P ≈ 125 MPa has a remarkably similar density and average structure as MDA, implying that MDA is not unique. Using the potential energy landscape formalism, we provide an intuitive qualitative understanding of the nature of LDA, HDA, and the IA generated at different pressures. In this view, LDA and HDA occupy specific and well-separated regions of the PEL; the IA prepared at P = 125 MPa is located in the intermediate region of the PEL that separates LDA and HDA.

Influence of surface nanostructure-induced innermost ion structuring on capacitance of carbon/ionic liquid double layers

Ionic liquids have drawn great interest as electrolytes for energy storage applications in which they form characteristic electrical double layers at electrode interfaces. For ionic liquids at carbon electrode interfaces, their double layers are subject to nanoscale structuring of the electrode surface, involving altered ion structure and interactions that significantly influence the double layer capacitance. In this regard, we investigate the modulation of ionic liquid double layers by electrode surface roughness and the resulting effects on the ion structure, interaction, and capacitance. We performed fixed voltage molecular dynamics simulations to compute the differential capacitance profiles for the ionic liquids [BMIm+][TFSI-] and [BMIm+][FSI-] at model carbon electrode interfaces with the surface channel width at subnanometer and nanometer scales. We find that both [BMIm+][TFSI-] and [BMIm+][FSI-] exhibit enhanced differential capacitance for the electrode surface with a subnanometer channel width relative to the flat graphene surface, but the most pronounced enhancements for these two ionic liquids unexpectedly appear at different applied potential regimes. For [BMIm+][TFSI-], the nanostructured electrode shows significant enhancement of capacitance at high positive potential. For [BMIm+][FSI-], on the other hand, this enhancement is small at positive polarization but noticeable at low negative potential. We demonstrate that differences in these capacitance trends is due to differences in ion correlation that arise from a steric constraint of nanostructured electrode surface on the voltage-mediated restructuring of ions closest to the electrode interface. For example, the TFSI- and FSI- anions tend to structure with their charged and nonpolar groups in contact with the positive electrode surface when the constraint on these close-contact anions is relaxed. This anion structuring largely retains the cation association near the nanostructured electrode, resulting in only a slight increase in capacitance at positive polarization. Our simulations highlight the sensitive dependence of the innermost ion structure on the electrode surface nanostructure and applied voltage and the resulting influence on ion correlation and capacitance of ionic liquid double layers.

Green-engineered clay- and carbon-based composite materials for the adsorption of benzene from air

Benzene is a carcinogenic volatile organic compound (VOC) that is ubiquitously detected in enclosed spaces due to emissions from cooking activities, building materials, and cleaning products. To remove benzene and other VOCs from indoor air and protect public health, traditional fabric filters have been modified to contain activated carbons to enhance the filtration efficacy. In this study, composites derived from natural clay minerals and activated carbon were individually green-engineered with chlorophylls and were attached to the surface of filter materials. These systems were assessed for their adsorption of benzene from air using in vitro and in silico methods. Isothermal, thermodynamic, and kinetic experiments indicated that all green-engineered composites had improved binding profiles for benzene, as demonstrated by increased binding affinities (Kf ≥ 900 vs 472) and lower values of Gibbs free energy (ΔG = -16.8 vs -15.2) compared to activated carbon. Adsorption of benzene to all composites was achieved quickly (< 30 min), and the green-engineered composites also showed low levels of desorption (≤ 25%). While free chlorophyll is known to be photosensitive, chlorophylls in the green-engineered composites showed photostability and maintained high binding rates (≥ 70%). Additionally, the in silico simulations demonstrated the significant contribution of chlorophyll for the overall binding of benzene in clay systems and that chlorophyll could contribute to benzene binding in the carbon-based systems. Together, these studies indicated that novel, green-engineered composite materials can be effective filter sorbents to enhance the removal of benzene from air.

Adsorption and removal of polystyrene nanoplastics from water by green-engineered clays

Human exposure to micro- and nanoplastics (MNPs) commonly occurs through the consumption of contaminated drinking water. Among these, polystyrene (PS) is well-characterized and is one of the most abundant MNPs, accounting for 10 % of total plastics. Previous studies have focused on carbonaceous materials to remove MNPs by filtration, but most of the work has involved microplastics since nanoplastics (NPs) are smaller in size and more difficult to measure and remove. To address this need, green-engineered chlorophyll-amended sodium and calcium montmorillonites (SMCH and CMCH) were tested for their ability to bind and detoxify parent and fluorescently labeled PSNP using in vitro, in silico, and in vivo assays. In vitro dosimetry, isothermal analyses, thermodynamics, and adsorption/desorption kinetic models demonstrated 1) high binding capacities (173-190 g/kg), 2) high affinities (103), and 3) chemisorption as suggested by low desorption (≤42 %) and high Gibbs free energy and enthalpy (>|-20| kJ/mol) in the Langmuir and pseudo-second-order models. Computational dynamics simulations for 30 and 40 monomeric units of PSNP depicted that chlorophyll amendments increased the binding percentage and contributed to the sustained binding. Also, 64 % of PSNP bind to both the head and tail of chlorophyll aggregates, rather than the head or tail only. Fluorescent PSNP at 100 nm and 30 nm that were exposed to Hydra vulgaris showed concentration-dependent toxicity at 20-100 µg/mL. Importantly, the inclusion of 0.05-0.3 % CMCH and SMCH significantly (p ≤ 0.01) and dose-dependently reduced PSNP toxicity in morphological changes and feeding rate. The bioassay validated the in vitro and in silico predictions about adsorption efficacy and mechanisms and suggested that CMCH and SMCH are efficacious binders for PSNP in water.

Toward Accurate Simulation of Coupling between Protein Secondary Structure and Phase Separation

Intrinsically disordered proteins (IDPs) frequently mediate phase separation that underlies the formation of a biomolecular condensate. Together with theory and experiment, efficient coarse-grained (CG) simulations have been instrumental in understanding the sequence-specific phase separation of IDPs. However, the widely used Cα-only models are limited in capturing the peptide nature of IDPs, particularly backbone-mediated interactions and effects of secondary structures, in phase separation. Here, we describe a hybrid resolution (HyRes) protein model toward a more accurate description of the backbone and transient secondary structures in phase separation. With an atomistic backbone and coarse-grained side chains, HyRes can semiquantitatively capture the residue helical propensity and overall chain dimension of monomeric IDPs. Using GY-23 as a model system, we show that HyRes is efficient enough for the direct simulation of spontaneous phase separation and, at the same time, appears accurate enough to resolve the effects of single His to Lys mutations. HyRes simulations also successfully predict increased β-structure formation in the condensate, consistent with available experimental CD data. We further utilize HyRes to study the phase separation of TPD-43, where several disease-related mutants in the conserved region (CR) have been shown to affect residual helicities and modulate the phase separation propensity as measured by the saturation concentration. The simulations successfully recapitulate the effect of these mutants on the helicity and phase separation propensity of TDP-43 CR. Analyses reveal that the balance between backbone and side chain-mediated interactions, but not helicity itself, actually determines phase separation propensity. These results support that HyRes represents an effective protein model for molecular simulation of IDP phase separation and will help to elucidate the coupling between transient secondary structures and phase separation.

Computational Characterization of the Reactivity of Compound I in Unspecific Peroxygenases

Unspecific peroxygenases (UPOs) are emerging as promising biocatalysts for selective oxyfunctionalization of unactivated C-H bonds. However, their potential in large-scale synthesis is currently constrained by suboptimal chemical selectivity. Improving the selectivity of UPOs requires a deep understanding of the molecular basis of their catalysis. Recent molecular simulations have sought to unravel UPO's selectivity and inform their design principles. However, most of these studies focused on substrate-binding poses. Few researchers have investigated how the reactivity of CpdI, the principal oxidizing intermediate in the catalytic cycle, influences selectivity in a realistic protein environment. Moreover, the influence of protein electrostatics on the reaction kinetics of CpdI has also been largely overlooked. To bridge this gap, we used multiscale simulations to interpret the regio- and enantioselective hydroxylation of the n-heptane substrate catalyzed by Agrocybe aegerita UPO (AaeUPO). We comprehensively characterized the energetics and kinetics of the hydrogen atom-transfer (HAT) step, initiated by CpdI, and the subsequent oxygen rebound step forming the product. Notably, our approach involved both free energy and potential energy evaluations in a quantum mechanics/molecular mechanics (QM/MM) setting, mitigating the dependence of results on the choice of initial conditions. These calculations illuminate the thermodynamics and kinetics of the HAT and oxygen rebound steps. Our findings highlight that both the conformational selection and the distinct chemical reactivity of different substrate hydrogen atoms together dictate the regio- and enantio-selectivity. Building on our previous study of CpdI's formation in AaeUPO, our results indicate that the HAT step is the rate-limiting step in the overall catalytic cycle. The subsequent oxygen rebound step is swift and retains the selectivity determined by the HAT step. We also pinpointed several polar and charged amino acid residues whose electrostatic potentials considerably influence the reaction barrier of the HAT step. Notably, the Glu196 residue is pivotal for both the CpdI's formation and participation in the HAT step. Our research offers in-depth insights into the catalytic cycle of AaeUPO, which will be instrumental in the rational design of UPOs with enhanced properties.

Backmapping with Mapping and Isomeric Information

I present a powerful and flexible backmapping tool named Multiscale Simulation Tool (mstool) that converts a coarse-grained (CG) system into all-atom (AA) resolution and only requires AA to CG mapping and isomeric information (cis/trans/dihedral/chiral). The backmapping procedure includes two simple steps: (a) AA atoms are randomly placed near the corresponding CG beads according to the provided mapping scheme. (b) Energy minimization is performed with two modifications in the AA force field (FF). First, nonbonded interactions are replaced with cosine functions to ensure the numerical stability. Second, additional torsions are imposed to maintain the molecules' isomeric properties. To test the simplicity and robustness of the tool, I backmapped multiple membrane and protein CG structures into AA resolution, including a four-bead CG lipid model (resolution increased by a factor of 34) without using intermediate resolution. The tool is freely available at github.com/ksy141/mstool.

Interactive molecular dynamics in virtual reality for modelling materials and catalysts

Interactive molecular dynamics simulation in virtual reality (iMD-VR) is emerging as a promising technique in molecular science. Here, we demonstrate its use in a range of fifteen applications in materials science and heterogeneous catalysis. In this work, the iMD-VR package Narupa is used with the MD package, DL_POLY [1]. We show how iMD-VR can be used to: (i) investigate the mechanism of lithium fast ion conduction by directing the formation of defects showing that vacancy transport is favoured over interstitialcy mechanisms, and (ii) guide a molecule through a zeolite pore to explore diffusion within zeolites, examining in detail the motion of methyl n-hexanoate in H-ZSM-5 zeolite and identifying bottlenecks restricting diffusion. iMD-VR allows users to manipulate these systems intuitively, to drive changes in them and observe the resulting changes in structure and dynamics. We make these simulations available, as a resource for both teaching and research. All simulation files, with videos, can be found online (https://doi.org/10.5281/zenodo.8252314) and are provided as open-source material.

Effect of Surface Chemistry on the Electrical Double Layer in a Long-Chain Ionic Liquid

Room temperature ionic liquids (ILs) can create a strong accumulation of charges at solid interfaces by forming a very thin and dense electrical double layer (EDL). The structure of this EDL has important consequences in numerous applications involving ILs, for example, in supercapacitors, sensors, and lubricants, by impacting the interfacial capacitance, the charge carrier density of semiconductors, as well as the frictional properties of the interfaces. We have studied the interfacial structure of a long chain imidazolium-based IL (1-octyl-3-methylimidazolium dicyanamide) on several substrates: mica, silica, silicon, and molybdenum disulfide (MoS2), using atomic force microscopy (AFM) experiments and molecular dynamics (MD) simulations. We have observed 3 types of interfacial structures for the same IL, depending on the chemistry of the substrate and the water content, showing that the EDL structure is not an intrinsic property of the IL. We evidenced that at a low water content, neutral and apolar (thus hydrophobic) substrates promote a thin layer structure, where the ions are oriented parallel to the substrate and cations and anions are mixed in each layer. In contrast, a strongly charged (thus hydrophilic) substrate yields an extended structuration into several bilayers, while a heterogeneous layering with loose bilayer regions was observed on an intermediate polar and weakly charged substrate and on an apolar one at a high bulk water content. In the latter case, water contamination favors the formation of bilayer patches by promoting the segregation of the long chain IL into polar and apolar domains.

Mobile barrier mechanisms for Na+-coupled symport in an MFS sugar transporter

While many 3D structures of cation-coupled transporters have been determined, the mechanistic details governing the obligatory coupling and functional regulations still remain elusive. The bacterial melibiose transporter (MelB) is a prototype of the Na+-coupled major facilitator superfamily transporters. With a conformational nanobody (Nb), we determined a low-sugar affinity inward-facing Na+-bound cryoEM structure. Collectively with the available outward-facing sugar-bound structures, both the outer and inner barriers were localized. The N- and C-terminal residues of the inner barrier contribute to the sugar selectivity pocket. When the inner barrier is broken as shown in the inward-open conformation, the sugar selectivity pocket is also broken. The binding assays by isothermal titration calorimetry revealed that this inward-facing conformation trapped by the conformation-selective Nb exhibited a greatly decreased sugar-binding affinity, suggesting the mechanisms for the substrate intracellular release and accumulation. While the inner/outer barrier shift directly regulates the sugar-binding affinity, it has little or no effect on the cation binding, which is also supported by molecular dynamics simulations. Furthermore, the use of this Nb in combination with the hydron/deuterium exchange mass spectrometry allowed us to identify dynamic regions; some regions are involved in the functionally important inner barrier-specific salt-bridge network, which indicates their critical roles in the barrier switching mechanisms for transport. These complementary results provided structural and dynamic insights into the mobile barrier mechanism for cation-coupled symport.

Binding of small molecule inhibitors to RNA polymerase-Spt5 complex impacts RNA and DNA stability

Spt5 is an elongation factor that associates with RNA polymerase II (Pol II) during transcription and has important functions in promoter-proximal pausing and elongation processivity. Spt5 was also recognized for its roles in the transcription of expanded-repeat genes that are related to neurodegenerative diseases. Recently, a set of Spt5-Pol II small molecule inhibitors (SPIs) were reported, which selectively inhibit mutant huntingtin gene transcription. Inhibition mechanisms as well as interaction sites of these SPIs with Pol II and Spt5 are not entirely known. In this study, we predicted the binding sites of three selected SPIs at the Pol II-Spt5 interface by docking and molecular dynamics simulations. Two molecules out of three demonstrated strong binding with Spt5 and Pol II, while the other molecule was more loosely bound and sampled multiple binding sites. Strongly bound SPIs indirectly affected RNA and DNA dynamics at the exit site as DNA became more flexible while RNA was stabilized by increased interactions with Spt5. Our results suggest that the transcription inhibition mechanism induced by SPIs can be related to Spt5-nucleic acid interactions, which were altered to some extent with strong binding of SPIs.

Alchemical Free Energy Calculations on Membrane-Associated Proteins

Membrane proteins have diverse functions within cells and are well-established drug targets. The advances in membrane protein structural biology have revealed drug and lipid binding sites on membrane proteins, while computational methods such as molecular simulations can resolve the thermodynamic basis of these interactions. Particularly, alchemical free energy calculations have shown promise in the calculation of reliable and reproducible binding free energies of protein-ligand and protein-lipid complexes in membrane-associated systems. In this review, we present an overview of representative alchemical free energy studies on G-protein-coupled receptors, ion channels, transporters as well as protein-lipid interactions, with emphasis on best practices and critical aspects of running these simulations. Additionally, we analyze challenges and successes when running alchemical free energy calculations on membrane-associated proteins. Finally, we highlight the value of alchemical free energy calculations calculations in drug discovery and their applicability in the pharmaceutical industry.

Accessible Molecular System Creator: Building Molecular Configurations Based on the Inaccessible Molecular Volume and Accessible Molecular Surface via Static Monte Carlo Sampling

Monte Carlo (MC) stochastic sampling is a powerful tool in classical molecular simulations that directly connects the observable macroscopic properties of matter and the underlying atomistic interactions. This connection operates within the framework of the statistical mechanics proposed by Gibbs. Most MC simulations are "dynamic," creating statistical ensembles of microstates via a Markovian chain, where each microstate in the ensemble depends only on its previous microstate. Herein, we re-examine an alternative form of MC that generates ensemble members through a "static" approach, building molecular systems stepwise. The basic theory for such an approach traces back to Rosenbluth and Rosenbluth, who proposed "static" stepwise sampling of a polymeric chain. It is almost as old as the Metropolis importance sampling approach used in dynamic MC, although the latter has been considerably more popular than the former. Herein, we address the main obstacle in static MC that has hindered the widespread adoption of Rosenbluth-based approaches in atomistic simulations. The obstacle lies in mapping the molecular accessible volume for adding a molecule in a Rosenbluth-like static sampling of atomistic configurations. We demonstrate a breakthrough by leveraging the ability to analytically map the inaccessible molecular volume and the accessible molecular surface owing to interatomically excluded volume interactions. This advance substantially enhances the ability to create molecular samples using a Rosenbluth-like static building process. The proposed approach can be used as a tool for creating initial configurations in MC or molecular dynamics simulations─a field where Rosenbluth-like static building has been applied. Additionally, this approach can be used as the first step in a perturbation scheme that accurately estimates free energy differences by estimating the chemical work related to molecule addition, removal, or reinsertion within the context of free energy perturbation schemes employed in molecular simulations.

Complex Conformational Space of the RNA Polymerase II C-Terminal Domain upon Phosphorylation

Intrinsically disordered proteins (IDPs) have been closely studied during the past decade due to their importance in many biological processes. The disordered nature of this group of proteins makes it difficult to observe its full span of the conformational space using either experimental or computational studies. In this article, we explored the conformational space of the C-terminal domain (CTD) of RNA polymerase II (Pol II), which is also an intrinsically disordered low complexity domain, using enhanced sampling methods. We provided a detailed conformational analysis of model systems of CTD with different lengths; first with the last 44 residues of the human CTD sequence and finally the CTD model with 2-heptapeptide repeating units. We then investigated the effects of phosphorylation on CTD conformations by performing simulations at different phosphorylated states. We obtained broad conformational spaces in nonphosphorylated CTD models, and phosphorylation has complex effects on the conformations of the CTD. These complex effects depend on the length of the CTD, spacing between the multiple phosphorylation sites, ion coordination, and interactions with the nearby residues.

Granger causality based on vector time series and quaternion algebra with possible applications to molecular dynamics data analysis

Causal analysis plays a significant role in physics, chemistry, and biology. Dynamics of complex (bio)molecular and nanosystems, from the microscopic to the macroscopic scale, are characterized by time-dependent vectors such as positions, forces, momenta, angular momenta, or torques. Identification and analysis of causal relationships between these time-dependent signals is an important problem in the multidimensional time-series analysis and is of great practical importance in describing the properties of such dynamical systems, and to understanding their functionality. For linear stochastic systems characterized by multidimensional scalar signals, Granger proposed a simple procedure to detect causal relationships, called Granger causality. In this study we extended this formalism to vector signals representing physical vector quantities. For this purpose, we used quaternion algebra, where vector signals are treated as time-dependent quaternions. The developed analytical model is based on the autoregressive formalism. This formalism (Q-MVAR) and its numerical implementation were validated using two simple dynamic models: a rigid body model represented by a benzenelike molecular fragment, interacting with a short-range harmonic potential with a wall, as well as a system of three model atomic balls moving inside a soft spherical surface and interacting with long range electrostatic forces. Although the motivation to these studies was the analysis of classical motions in complex (bio)molecular systems, described with a mechanical model and based on molecular dynamics (MD) simulations, in particular coarse-grained ones, it should be noted that the developed extended formalism can be applied to any system composed of many rigid elements that interact with arbitrary potentials and are characterized by complex internal motions. A description of the detailed procedure for calculating causality measures is provided in the Appendices of the Supplemental Material. This formalism and the prototype of its numerical implementation can be further developed and applied in many different fields of physical, natural, and engineering sciences.

Understanding the Multifaceted Mechanism of Compound I Formation in Unspecific Peroxygenases through Multiscale Simulations

Unspecific peroxygenases (UPOs) can selectively oxyfunctionalize unactivated hydrocarbons by using peroxides under mild conditions. They circumvent the oxygen dilemma faced by cytochrome P450s and exhibit greater stability than the latter. As such, they hold great potential for industrial applications. A thorough understanding of their catalysis is needed to improve their catalytic performance. However, it remains elusive how UPOs effectively convert peroxide to Compound I (CpdI), the principal oxidizing intermediate in the catalytic cycle. Previous computational studies of this process primarily focused on heme peroxidases and P450s, which have significant differences in the active site from UPOs. Additionally, the roles of peroxide unbinding in the kinetics of CpdI formation, which is essential for interpreting existing experiments, have been understudied. Moreover, there has been a lack of free energy characterizations with explicit sampling of protein and hydration dynamics, which is critical for understanding the thermodynamics of the proton transport (PT) events involved in CpdI formation. To bridge these gaps, we employed multiscale simulations to comprehensively characterize the CpdI formation in wild-type UPO from Agrocybe aegerita (AaeUPO). Extensive free energy and potential energy calculations were performed in a quantum mechanics/molecular mechanics setting. Our results indicate that substrate-binding dehydrates the active site, impeding the PT from H2O2 to a nearby catalytic base (Glu196). Furthermore, the PT is coupled with considerable hydrogen bond network rearrangements near the active site, facilitating subsequent O-O bond cleavage. Finally, large unbinding free energy barriers kinetically stabilize H2O2 at the active site. These findings reveal a delicate balance among PT, hydration dynamics, hydrogen bond rearrangement, and cosubstrate unbinding, which collectively enable efficient CpdI formation. Our simulation results are consistent with kinetic measurements and offer new insights into the CpdI formation mechanism at atomic-level details, which can potentially aid the design of next-generation biocatalysts for sustainable chemical transformations of feedstocks.

Computing Free Energies of Fold-Switching Proteins Using MELD x MD

Some proteins are conformational switches, able to transition between relatively different conformations. To understand what drives them requires computing the free-energy difference ΔGAB between their stable states, A and B. Molecular dynamics (MD) simulations alone are often slow because they require a reaction coordinate and must sample many transitions in between. Here, we show that modeling employing limited data (MELD) x MD on known endstates A and B is accurate and efficient because it does not require passing over barriers or knowing reaction coordinates. We validate this method on two problems: (1) it gives correct relative populations of α and β conformers for small designed chameleon sequences of protein G; and (2) it correctly predicts the conformations of the C-terminal domain (CTD) of RfaH. Free-energy methods like MELD x MD can often resolve structures that confuse machine-learning (ML) methods.

G-Quadruplex-Filtered Selective Ion-to-Ion Current Amplification for Non-Invasive Ion Monitoring in Real Time

Living cells efflux intracellular ions for maintaining cellular life, so intravital measurements of specific ion signals are of significant importance for studying cellular functions and pharmacokinetics. In this work, de novo synthesis of artificial K+ -selective membrane and its integration with polyelectrolyte hydrogel-based open-junction ionic diode (OJID) is demonstrated, achieving a real-time K+ -selective ion-to-ion current amplification in complex bioenvironments. By mimicking biological K+ channels and nerve impulse transmitters, in-line K+ -binding G-quartets are introduced across freestanding lipid bilayers by G-specific hexylation of monolithic G-quadruplex, and the pre-filtered K+ flow is directly converted to amplified ionic currents by the OJID with a fast response time at 100 ms intervals. By the synergistic combination of charge repulsion, sieving, and ion recognition, the synthetic membrane allows K+ transport exclusively without water leakage; it is 250× and 17× more permeable toward K+ than monovalent anion, Cl- , and polyatomic cation, N-methyl-d-glucamine+ , respectively. The molecular recognition-mediated ion channeling provides a 500% larger signal for K+ as compared to Li+ (0.6× smaller than K+ ) despite the same valence. Using the miniaturized device, non-invasive, direct, and real-time K+ efflux monitoring from living cell spheroids is achieved with minimal crosstalk, specifically in identifying osmotic shock-induced necrosis and drug-antidote dynamics.

Structural basis of transcription: RNA Polymerase II substrate binding and metal coordination at 3.0 Å using a free-electron laser

Catalysis and translocation of multi-subunit DNA-directed RNA polymerases underlie all cellular mRNA synthesis. RNA polymerase II (Pol II) synthesizes eukaryotic pre-mRNAs from a DNA template strand buried in its active site. Structural details of catalysis at near atomic resolution and precise arrangement of key active site components have been elusive. Here we present the free electron laser (FEL) structure of a matched ATP-bound Pol II, revealing the full active site interaction network at the highest resolution to date, including the trigger loop (TL) in the closed conformation, bonafide occupancy of both site A and B Mg2+, and a putative third (site C) Mg2+ analogous to that described for some DNA polymerases but not observed previously for cellular RNA polymerases. Molecular dynamics (MD) simulations of the structure indicate that the third Mg2+ is coordinated and stabilized at its observed position. TL residues provide half of the substrate binding pocket while multiple TL/bridge helix (BH) interactions induce conformational changes that could propel translocation upon substrate hydrolysis. Consistent with TL/BH communication, a FEL structure and MD simulations of the hyperactive Rpb1 T834P bridge helix mutant reveals rearrangement of some active site interactions supporting potential plasticity in active site function and long-distance effects on both the width of the central channel and TL conformation, likely underlying its increased elongation rate at the expense of fidelity.

Incorporating physics to overcome data scarcity in predictive modeling of protein function: A case study of BK channels

Machine learning has played transformative roles in numerous chemical and biophysical problems such as protein folding where large amount of data exists. Nonetheless, many important problems remain challenging for data-driven machine learning approaches due to the limitation of data scarcity. One approach to overcome data scarcity is to incorporate physical principles such as through molecular modeling and simulation. Here, we focus on the big potassium (BK) channels that play important roles in cardiovascular and neural systems. Many mutants of BK channel are associated with various neurological and cardiovascular diseases, but the molecular effects are unknown. The voltage gating properties of BK channels have been characterized for 473 site-specific mutations experimentally over the last three decades; yet, these functional data by themselves remain far too sparse to derive a predictive model of BK channel voltage gating. Using physics-based modeling, we quantify the energetic effects of all single mutations on both open and closed states of the channel. Together with dynamic properties derived from atomistic simulations, these physical descriptors allow the training of random forest models that could reproduce unseen experimentally measured shifts in gating voltage, ∆V1/2, with a RMSE ~ 32 mV and correlation coefficient of R ~ 0.7. Importantly, the model appears capable of uncovering nontrivial physical principles underlying the gating of the channel, including a central role of hydrophobic gating. The model was further evaluated using four novel mutations of L235 and V236 on the S5 helix, mutations of which are predicted to have opposing effects on V1/2 and suggest a key role of S5 in mediating voltage sensor-pore coupling. The measured ∆V1/2 agree quantitatively with prediction for all four mutations, with a high correlation of R = 0.92 and RMSE = 18 mV. Therefore, the model can capture nontrivial voltage gating properties in regions where few mutations are known. The success of predictive modeling of BK voltage gating demonstrates the potential of combining physics and statistical learning for overcoming data scarcity in nontrivial protein function prediction.

Accurate Simulation of Coupling between Protein Secondary Structure and Liquid-Liquid Phase Separation

Intrinsically disordered proteins (IDPs) frequently mediate liquid-liquid phase separation (LLPS) that underlies the formation of membraneless organelles. Together with theory and experiment, efficient coarse-grained (CG) simulations have been instrumental in understanding sequence-specific phase separation of IDPs. However, the widely-used Cα-only models are severely limited in capturing the peptide nature of IDPs, including backbone-mediated interactions and effects of secondary structures, in LLPS. Here, we describe a hybrid resolution (HyRes) protein model for accurate description of the backbone and transient secondary structures in LLPS. With an atomistic backbone and coarse-grained side chains, HyRes accurately predicts the residue helical propensity and chain dimension of monomeric IDPs. Using GY-23 as a model system, we show that HyRes is efficient enough for direct simulation of spontaneous phase separation, and at the same time accurate enough to resolve the effects of single mutations. HyRes simulations also successfully predict increased beta-sheet formation in the condensate, consistent with available experimental data. We further utilize HyRes to study the phase separation of TPD-43, where several disease-related mutants in the conserved region (CR) have been shown to affect residual helicities and modulate LLPS propensity. The simulations successfully recapitulate the effect of these mutants on the helicity and LLPS propensity of TDP-43 CR. Analyses reveal that the balance between backbone and sidechain-mediated interactions, but not helicity itself, actually determines LLPS propensity. We believe that the HyRes model represents an important advance in the molecular simulation of LLPS and will help elucidate the coupling between IDP transient secondary structures and phase separation.

Aggregation-Induced Structural Symmetry Breaking Promotes Charge Separation for Efficient Photocatalytic Hydrogen Production

Recently, organic semiconductors have received much attention in the field of photocatalysis due to their tunable physicochemical properties. However, organic semiconductor photocatalysts typically suffer from severe charge recombination due to high exciton binding energy. Herein, we found that aggregation of pyrene results in a red-shift of the light absorption from UV to visible light region. Importantly, the aggregation can induce dipole polarization by spontaneous structural symmetry breaking, thus significantly accelerating the separation and transfer of charge carriers. As a result, the pyrene aggregates display enhanced hydrogen photosynthesis activity. Furthermore, the noncovalent interactions allow rational design of physicochemical and electronic properties of pyrene aggregates, further strengthening the charge separation and photocatalytic activity of aggregates. The quantum yield of pyrene aggregates for hydrogen production highly reaches 20.77 % at 400 nm. Moreover, we have also observed pyrene analogues (1-hydroxypyrene, 1-nitropyrene and perylene) after aggregation all display large dipole moments induced by structural symmetry breaking and therefore accelerate the separation of charge carriers, confirming its general principle. This work highlights the achievement of using aggregation-induced structural symmetry breaking to enable the separation and transfer of charge carriers.

Data scheme and data format for transferable force fields for molecular simulation

A generalized data scheme for transferable classical force fields used in molecular simulations, i.e. molecular dynamics and Monte Carlo simulation, is presented. The data scheme is implemented in an SQL-based data format. The data scheme and data format is machine readable, re-usable, and interoperable. A transferable force field is a chemical construction plan specifying intermolecular and intramolecular interactions between different types of atoms or different chemical groups and can be used for building a model for a given component. The data scheme proposed in this work (named TUK-FFDat) formalizes digitally these chemical construction plans, i.e. transferable force fields. It can be applied to all-atom as well as united-atom transferable force fields. The general applicability of the data scheme is demonstrated for different types of force fields (TraPPE, OPLS-AA, and Potoff). Furthermore, conversion tools for translating the data scheme between .xls spread sheet format and the SQL-based data format are provided. The data format can readily be integrated in existing workflows, simulation engines, and force field databases as well as for linking such.

EnGens: a computational framework for generation and analysis of representative protein conformational ensembles

Proteins are dynamic macromolecules that perform vital functions in cells. A protein structure determines its function, but this structure is not static, as proteins change their conformation to achieve various functions. Understanding the conformational landscapes of proteins is essential to understand their mechanism of action. Sets of carefully chosen conformations can summarize such complex landscapes and provide better insights into protein function than single conformations. We refer to these sets as representative conformational ensembles. Recent advances in computational methods have led to an increase in the number of available structural datasets spanning conformational landscapes. However, extracting representative conformational ensembles from such datasets is not an easy task and many methods have been developed to tackle it. Our new approach, EnGens (short for ensemble generation), collects these methods into a unified framework for generating and analyzing representative protein conformational ensembles. In this work, we: (1) provide an overview of existing methods and tools for representative protein structural ensemble generation and analysis; (2) unify existing approaches in an open-source Python package, and a portable Docker image, providing interactive visualizations within a Jupyter Notebook pipeline; (3) test our pipeline on a few canonical examples from the literature. Representative ensembles produced by EnGens can be used for many downstream tasks such as protein-ligand ensemble docking, Markov state modeling of protein dynamics and analysis of the effect of single-point mutations.

pyCHARMM: Embedding CHARMM Functionality in a Python Framework

CHARMM is rich in methodology and functionality as one of the first programs addressing problems of molecular dynamics and modeling of biological macromolecules and their partners, e.g., small molecule ligands. When combined with the highly developed CHARMM parameters for proteins, nucleic acids, small molecules, lipids, sugars, and other biologically relevant building blocks, and the versatile CHARMM scripting language, CHARMM has been a trendsetting platform for modeling studies of biological macromolecules. To further enhance the utility of accessing and using CHARMM functionality in increasingly complex workflows associated with modeling biological systems, we introduce pyCHARMM, Python bindings, functions, and modules to complement and extend the extensive set of modeling tools and methods already available in CHARMM. These include access to CHARMM function-generated variables associated with the system (psf), coordinates, velocities and forces, atom selection variables, and force field related parameters. The ability to augment CHARMM forces and energies with energy terms or methods derived from machine learning or other sources, written in Python, CUDA, or OpenCL and expressed as Python callable routines is introduced together with analogous functions callable during dynamics calculations. Integration of Python-based graphical engines for visualization of simulation models and results is also accessible. Loosely coupled parallelism is available for workflows such as free energy calculations, using MBAR/TI approaches or high-throughput multisite λ-dynamics (MSλD) free energy methods, string path optimization calculations, replica exchange, and molecular docking with a new Python-based CDOCKER module. CHARMM accelerated platform kernels through the CHARMM/OpenMM API, CHARMM/DOMDEC, and CHARMM/BLaDE API are also readily integrated into this Python framework. We anticipate that pyCHARMM will be a robust platform for the development of comprehensive and complex workflows utilizing Python and its extensive functionality as well as an optimal platform for users to learn molecular modeling methods and practices within a Python-friendly environment such as Jupyter Notebooks.

Incorporating physics to overcome data scarcity in predictive modeling of protein function: a case study of BK channels

Machine learning has played transformative roles in numerous chemical and biophysical problems such as protein folding where large amount of data exists. Nonetheless, many important problems remain challenging for data-driven machine learning approaches due to the limitation of data scarcity. One approach to overcome data scarcity is to incorporate physical principles such as through molecular modeling and simulation. Here, we focus on the big potassium (BK) channels that play important roles in cardiovascular and neural systems. Many mutants of BK channel are associated with various neurological and cardiovascular diseases, but the molecular effects are unknown. The voltage gating properties of BK channels have been characterized for 473 site-specific mutations experimentally over the last three decades; yet, these functional data by themselves remain far too sparse to derive a predictive model of BK channel voltage gating. Using physics-based modeling, we quantify the energetic effects of all single mutations on both open and closed states of the channel. Together with dynamic properties derived from atomistic simulations, these physical descriptors allow the training of random forest models that could reproduce unseen experimentally measured shifts in gating voltage, ΔV 1/2 , with a RMSE ∼ 32 mV and correlation coefficient of R ∼ 0.7. Importantly, the model appears capable of uncovering nontrivial physical principles underlying the gating of the channel, including a central role of hydrophobic gating. The model was further evaluated using four novel mutations of L235 and V236 on the S5 helix, mutations of which are predicted to have opposing effects on V 1/2 and suggest a key role of S5 in mediating voltage sensor-pore coupling. The measured ΔV 1/2 agree quantitatively with prediction for all four mutations, with a high correlation of R = 0.92 and RMSE = 18 mV. Therefore, the model can capture nontrivial voltage gating properties in regions where few mutations are known. The success of predictive modeling of BK voltage gating demonstrates the potential of combining physics and statistical learning for overcoming data scarcity in nontrivial protein function prediction.

Author Summary

Deep machine learning has brought many exciting breakthroughs in chemistry, physics and biology. These models require large amount of training data and struggle when the data is scarce. The latter is true for predictive modeling of the function of complex proteins such as ion channels, where only hundreds of mutational data may be available. Using the big potassium (BK) channel as a biologically important model system, we demonstrate that a reliable predictive model of its voltage gating property could be derived from only 473 mutational data by incorporating physics-derived features, which include dynamic properties from molecular dynamics simulations and energetic quantities from Rosetta mutation calculations. We show that the final random forest model captures key trends and hotspots in mutational effects of BK voltage gating, such as the important role of pore hydrophobicity. A particularly curious prediction is that mutations of two adjacent residues on the S5 helix would always have opposite effects on the gating voltage, which was confirmed by experimental characterization of four novel mutations. The current work demonstrates the importance and effectiveness of incorporating physics in predictive modeling of protein function with scarce data.

The Importance of Nuclear Quantum Effects on the Thermodynamic and Structural Properties of Low-Density Amorphous Ice: A Comparison with Hexagonal Ice

We study the nuclear quantum effects (NQE) on the thermodynamic properties of low-density amorphous ice (LDA) and hexagonal ice (I_h) at P = 0.1 MPa and T ≥ 25 K. Our results are based on path-integral molecular dynamics (PIMD) and classical MD simulations of H₂O and D₂O using the q-TIP4P/F water model. We show that the inclusion of NQE is necessary to reproduce the experimental properties of LDA and ice I_h. While MD simulations (no NQE) predict that the density ρ(T) of LDA and ice I_h increases monotonically upon cooling, PIMD simulations indicate the presence of a density maximum in LDA and ice I_h. MD and PIMD simulations also predict a qualitatively different T-dependence for the thermal expansion coefficient α_P(T) and bulk modulus B(T) of both LDA and ice I_h. Remarkably, the ρ(T), α_P(T), and B(T) of LDA are practically identical to those of ice I_h. The origin of the observed NQE is due to the delocalization of the H atoms, which is identical in LDA and ice I_h. H atoms delocalize considerably (over a distance ≈ 20-25% of the OH covalent-bond length) and anisotropically (preferentially perpendicular to the OH covalent bond), leading to less linear hydrogen bonds HB (larger HOO angles and longer OO separations) than observed in classical MD simulations.

Machine Learning Generation of Dynamic Protein Conformational Ensembles

Machine learning has achieved remarkable success across a broad range of scientific and engineering disciplines, particularly its use for predicting native protein structures from sequence information alone. However, biomolecules are inherently dynamic, and there is a pressing need for accurate predictions of dynamic structural ensembles across multiple functional levels. These problems range from the relatively well-defined task of predicting conformational dynamics around the native state of a protein, which traditional molecular dynamics (MD) simulations are particularly adept at handling, to generating large-scale conformational transitions connecting distinct functional states of structured proteins or numerous marginally stable states within the dynamic ensembles of intrinsically disordered proteins. Machine learning has been increasingly applied to learn low-dimensional representations of protein conformational spaces, which can then be used to drive additional MD sampling or directly generate novel conformations. These methods promise to greatly reduce the computational cost of generating dynamic protein ensembles, compared to traditional MD simulations. In this review, we examine recent progress in machine learning approaches towards generative modeling of dynamic protein ensembles and emphasize the crucial importance of integrating advances in machine learning, structural data, and physical principles to achieve these ambitious goals.

Cancer-cell-derived fumarate suppresses the anti-tumor capacity of CD8+ T cells in the tumor microenvironment

Metabolic alterations in the microenvironment significantly modulate tumor immunosensitivity, but the underlying mechanisms remain obscure. Here, we report that tumors depleted of fumarate hydratase (FH) exhibit inhibition of functional CD8⁺ T cell activation, expansion, and efficacy, with enhanced malignant proliferative capacity. Mechanistically, FH depletion in tumor cells accumulates fumarate in the tumor interstitial fluid, and increased fumarate can directly succinate ZAP70 at C96 and C102 and abrogate its activity in infiltrating CD8⁺ T cells, resulting in suppressed CD8⁺ T cell activation and anti-tumor immune responses in vitro and in vivo. Additionally, fumarate depletion by increasing FH expression strongly enhances the anti-tumor efficacy of anti-CD19 CAR T cells. Thus, these findings demonstrate a role for fumarate in controlling TCR signaling and suggest that fumarate accumulation in the tumor microenvironment (TME) is a metabolic barrier to CD8⁺ T cell anti-tumor function. And potentially, fumarate depletion could be an important strategy for tumor immunotherapy.

EnGens: a computational framework for generation and analysis of representative protein conformational ensembles

Proteins are dynamic macromolecules that perform vital functions in cells. A protein structure determines its function, but this structure is not static, as proteins change their conformation to achieve various functions. Understanding the conformational landscapes of proteins is essential to understand their mechanism of action. Sets of carefully chosen conformations can summarize such complex landscapes and provide better insights into protein function than single conformations. We refer to these sets as representative conformational ensembles. Recent advances in computational methods have led to an increase in number of available structural datasets spanning conformational landscapes. However, extracting representative conformational ensembles from such datasets is not an easy task and many methods have been developed to tackle it. Our new approach, EnGens (short for ensemble generation), collects these methods into a unified framework for generating and analyzing protein conformational ensembles. In this work we: (1) provide an overview of existing methods and tools for protein structural ensemble generation and analysis; (2) unify existing approaches in an open-source Python package, and a portable Docker image, providing interactive visualizations within a Jupyter Notebook pipeline; (3) test our pipeline on a few canonical examples found in the literature. Representative ensembles produced by EnGens can be used for many downstream tasks such as protein-ligand ensemble docking, Markov state modeling of protein dynamics and analysis of the effect of single-point mutations.

Improving drug-target affinity prediction via feature fusion and knowledge distillation

Rapid and accurate prediction of drug-target affinity can accelerate and improve the drug discovery process. Recent studies show that deep learning models may have the potential to provide fast and accurate drug-target affinity prediction. However, the existing deep learning models still have their own disadvantages that make it difficult to complete the task satisfactorily. Complex-based models rely heavily on the time-consuming docking process, and complex-free models lacks interpretability. In this study, we introduced a novel knowledge-distillation insights drug-target affinity prediction model with feature fusion inputs to make fast, accurate and explainable predictions. We benchmarked the model on public affinity prediction and virtual screening dataset. The results show that it outperformed previous state-of-the-art models and achieved comparable performance to previous complex-based models. Finally, we study the interpretability of this model through visualization and find it can provide meaningful explanations for pairwise interaction. We believe this model can further improve the drug-target affinity prediction for its higher accuracy and reliable interpretability.

Tools for analyzing protonation states and for tracing proton transfer pathways with examples from the Rb. sphaeroides photosynthetic reaction centers

Protons participate in many reactions. In proteins, protons need paths to move in and out of buried active sites. The vectorial movement of protons coupled to electron transfer reactions establishes the transmembrane electrochemical gradient used for many reactions, including ATP synthesis. Protons move through hydrogen bonded chains of waters and hydroxy side chains via the Grotthuss mechanism and by proton binding and release from acidic and basic residues. MCCE analysis shows that proteins exist in a large number of protonation states. Knowledge of the equilibrium ensemble can provide a rational basis for setting protonation states in simulations that fix them, such as molecular dynamics (MD). The proton path into the Q_B site in the bacterial reaction centers (RCs) of Rb. sphaeroides is analyzed by MD to provide an example of the benefits of using protonation states found by the MCCE program. A tangled web of side chains and waters link the cytoplasm to Q_B. MCCE analysis of snapshots from multiple trajectories shows that changing the input protonation state of a residue in MD biases the trajectory shifting the proton affinity of that residue. However, the proton affinity of some residues is more sensitive to the input structure. The proton transfer networks derived from different trajectories are quite robust. There are some changes in connectivity that are largely restricted to the specific residues whose protonation state is changed. Trajectories with Q_B^•- are compared with earlier results obtained with Q_B [Wei et. al Photosynthesis Research volume 152, pages153-165 (2022)] showing only modest changes. While introducing new methods the study highlights the difficulty of establishing the connections between protein conformation.

AmberMDrun: A Scripting Tool for Running Amber MD in an Easy Way

MD simulations have been widely applied and become a powerful tool in the field of biomacromolecule simulations and computer-aided drug design, etc., which can estimate binding free energy between receptor and ligand. However, the inputs and force field preparation for performing Amber MD is somewhat complicated, and challenging for beginners. To address this issue, we have developed a script for automatically preparing Amber MD input files, balancing the system, performing Amber MD for production, and predicting receptor-ligand binding free energy. This script is open-source, extensible and can support customization. The core code is written in C++ and has a Python interface, providing both efficient performance and convenience.

Nonadiabatic Forward Flux Sampling for Excited-State Rare Events

We present a rare event sampling scheme applicable to coupled electronic excited states. In particular, we extend the forward flux sampling (FFS) method for rare event sampling to a nonadiabatic version (NAFFS) that uses the trajectory surface hopping (TSH) method for nonadiabatic dynamics. NAFFS is applied to two dynamically relevant excited-state models that feature an avoided crossing and a conical intersection with tunable parameters. We investigate how nonadiabatic couplings, temperature, and reaction barriers affect transition rate constants in regimes that cannot be otherwise obtained with plain, traditional TSH. The comparison with reference brute-force TSH simulations for limiting cases of rareness shows that NAFFS can be several orders of magnitude cheaper than conventional TSH and thus represents a conceptually novel tool to extend excited-state dynamics to time scales that are able to capture rare nonadiabatic events.

Re-Balancing Replica Exchange with Solute Tempering for Sampling Dynamic Protein Conformations

Replica exchange with solute tempering (REST) is a highly effective variant of replica exchange for enhanced sampling in explicit solvent simulations of biomolecules. By scaling the Hamiltonian for a selected "solute" region of the system, REST effectively applies tempering only to the degrees of freedom of interest but not the rest of the system ("solvent"), allowing fewer replicas for covering the same temperature range. A key consideration of REST is how the solute-solvent interactions are scaled together with the solute-solute interactions. Here, we critically evaluate the performance of the latest REST2 protocol for sampling large-scale conformation fluctuations of intrinsically disordered proteins (IDPs). The results show that REST2 promotes artificial protein conformational collapse at high effective temperatures, which seems to be a designed feature originally to promote the sampling of reversible folding of small proteins. The collapse is particularly severe with larger IDPs, leading to replica segregation in the effective temperature space and hindering effective sampling of large-scale conformational changes. We propose that the scaling of the solute-solvent interactions can be treated as free parameters in REST, which can be tuned to control the solute conformational properties (e.g., chain expansion) at different effective temperatures and achieve more effective sampling. To this end, we derive a new REST3 protocol, where the strengths of the solute-solvent van der Waals interactions are recalibrated to reproduce the levels of protein chain expansion at high effective temperatures. The efficiency of REST3 is examined using two IDPs with nontrivial local and long-range structural features, including the p53 N-terminal domain and the kinase inducible transactivation domain of transcription factor CREB. The results suggest that REST3 leads to a much more efficient temperature random walk and improved sampling efficiency, which also further reduces the number of replicas required. Nonetheless, our analysis also reveals significant challenges of relying on tempering alone for sampling large-scale conformational fluctuations of disordered proteins. It is likely that more efficient sampling protocols will require incorporating more sophisticated Hamiltonian replica exchange schemes in addition to tempering.