The multipole and induced dipole (MPID) model provides a sophisticated description of electrostatics through permanent atomic multipoles up to the octopole level and point induced dipoles, but its high computational cost has limited its routine application to biomolecular simulations. Here, we present a GPU-accelerated implementation of the MPID model in OpenMM. The implementation produces energies and forces fully consistent with those obtained from CHARMM and enables routine microsecond-scale molecular dynamics simulations of explicitly solvated protein systems. Moreover, force field parameters were transferred directly from the Drude oscillator model to the MPID framework without reparameterization, yielding the MPID-2019 protein force field. Tests on five proteins show that MPID-2019 reproduces a range of NMR observables, including scalar couplings and relaxation order parameters. These results demonstrate that MPID-2019 is a valid polarizable force field readily applicable to protein simulations and that the effective equivalence of the induced dipole and Drude oscillator formalisms can be established at the biomolecular level.

Structures built for or generated from Molecular Dynamics simulations can often have problems such as wrong chiralities, bad parameter assignment, steric clashes, and/or distorted geometries. It is important to be able to detect such problems in a rapid and automated way since such issues can prevent subsequent simulations from running successfully. Here we detail a structure check procedure which focuses on checking for some of the more common structural issues, including atomic overlaps, distorted bonds, and ring intersections, as well as the strategies for parallelizing the check procedure so it can rapidly detect issues in structures with millions of atoms.

Enzymes catalyze complex chemical transformations with remarkable efficiency and selectivity, yet their atomistic mechanisms remain challenging to capture because conventional simulations trade accuracy for efficiency. Here we introduce a reactive machine learning/molecular mechanics (ML/MM) framework that bridges quantum chemistry with long-timescale sampling, enabling direct exploration of enzymatic transition states and free-energy landscapes. Coupled with metadynamics, this approach achieves nanosecond sampling of bond-forming reactions and quantitatively predicts activation barriers, mutational effects, and stereoselectivity. Applied to Diels-Alderases, the framework not only reproduces experimental activity and endo/exo preferences with sub-kcal mol-1 accuracy but also uncovers how pathway dynamics and local electrostatics preorganize substrates for selective outcomes. By uniting reactivity, conformational dynamics, and predictive power, this work establishes reactive ML/MM as a broadly applicable strategy for mechanistic enzymology and a foundation for the rational design of new biocatalysts.

Iron-sulfur (Fe-S) clusters are ubiquitous cofactors in metalloproteins, supporting essential biological functions such as electron transfer, enzymatic catalysis, and metabolic regulation. Despite their importance, large-scale identification of Fe-S proteins remains challenging due to limitations of experimental methods and inconsistent annotations in existing databases. To address this, we introduced FeSseqdb, a curated sequence-level database derived from the Protein Data Bank (PDB), in which Fe-S cluster-containing chains are systematically verified using atomic coordinates. By standardizing diverse ligand annotations, FeSseqdb provides a unified and reliable resource for Fe-S protein research. Building on this foundation, a machine learning framework was developed to predict Fe-S proteins using only sequence-derived features, including amino acid composition, cysteine-related metrics, and sequence length. Among the classifiers evaluated, the random forest model trained on data resampled with SVMSMOTE achieved the highest predictive performance, highlighting the discriminative power of simple sequence features. To elucidate the biological relevance of these features, explainable AI methods were applied to identify key sequence characteristics associated with Fe-S proteins. Cysteine frequency and spatial distribution, along with proline content, emerged as primary contributors, which is consistent with their known structural roles in cluster coordination. Additionally, serine, glutamic acid, and arginine were identified as secondary determinants, in line with their reported roles in redox and electrostatic environments surrounding metal cofactors. The inclusion of these biologically relevant features demonstrates the potential of sequence-based models not only for accurate prediction but also for uncovering functional insights that align with known biochemical principles. This approach provides a foundation for large-scale, sequence-based discovery of Fe-S proteins and supports future investigations into their functional diversity across proteomes.

Despite widespread assumptions of ZIF-8 stability in aqueous media, we demonstrate that the apparent endurance of ZIF-8 is largely misconceived as thermodynamic stability rather than correctly ascribed as kinetic stability under a certain experimental setup. In unbuffered solutions, the rapid release of zinc and 2-methylimidazole ligands causes an instantaneous pH spike, which inhibits further degradation. This self-protective mechanism is, however, nullified in buffered environments, where ZIF-8 undergoes a rapid and spontaneous decomposition. Our findings ascertain a critical thermodynamic instability of ZIF-8 under physiologically relevant conditions. The findings of this investigation necessitate a fundamental re-evaluation of ZIF-8 suitability for water-based applications, particularly drug delivery. Furthermore, the insights and methodology presented herein provide a crucial framework for predicting the stability of similar metal-organic frameworks in aqueous solutions under biologically relevant conditions.

Accurate prediction of half-maximal inhibitory concentration (IC50) values is critical for accelerating drug discovery, yet traditional quantitative structure-activity relationship (QSAR) models often have limited ability to capture both local structural patterns and global physicochemical properties essential for bioactivity. We developed a hybrid deep learning framework that integrates graph neural networks with explicit molecular descriptors to address this limitation. The model learns from molecular graphs encoding atomic and bond features while incorporating interpretable physicochemical properties and structural fingerprints. Trained and validated on 14,316 compounds across nine diverse biological targets including kinases, nuclear receptors, and proteases, our approach achieved an overall test R2 of 0.87, consistently outperforming previously reported methods by 6-42% across evaluated targets. The model demonstrated robust generalization with near-identical training and test performance, while maintaining partial interpretability through transparent descriptor contributions and attention mechanisms. By synergistically combining data-driven learning with domain knowledge, this hybrid framework offers improved accuracy and interpretability for structure-activity modeling, facilitating more efficient compound prioritization and optimization in early-stage drug discovery programs.

Glycosaminoglycans (GAGs) are long, anionic polysaccharides abundant in the extracellular matrix and lysosomes, where their electrostatic interactions with proteins are essential for biological function. Computational studies of GAG-containing systems remain challenging due to their significant charge density and conformational flexibility. Here we benchmark two widely used force-fields, ff14SB/GLYCAM06j-1 and CHARMM36m, for three experimentally characterized protein-GAG complexes. Both force fields reproduce the key structural features of protein-GAG interactions, while GAG dynamics depend on protein charge, with CHARMM36m favoring broader surface exploration for highly positively charged proteins and AMBER enhancing mobility for less charged systems. Although protein flexibility is similarly described, ff14SB/GLYCAM06j-1 samples a broader GAG conformational space, and dissociation free energy profiles diverge for highly anionic GAGs, but remain comparable for moderately sulfated systems. In addition, we performed molecular dynamics simulations for all systems using the ff14SB/GLYCAM06j-1, CHARMM36m, and ff19SB/GLYCAM06j-1 force fields in a 15 Å solvent box. Structural and energetic analyses revealed no significant impact of the solvent box size on the examined descriptors. These results establish practical benchmarks for accurate atomistic simulations of GAG-protein assemblies and will inform future developments in biomolecular force fields.

KRAS is a predominant oncogenic driver across multiple cancers and was long considered undruggable due to its high nucleotide affinity and lack of classical binding pockets. Although recent advances have led to covalent inhibitors such as Sotorasib and Adagrasib for the KRAS G12C mutation, effective therapies for other common variants, particularly KRAS G12D, which is highly prevalent in aggressive pancreatic cancers, remain limited. In this study, we employed machine learning approaches to identify potential inhibitors of KRAS G12D and G12C by screening FDA-approved compounds curated from the ChEMBL database. Random Forest and Neural Network models were trained on bioactivity data from three BindingDB data sets: wild-type KRAS GTPase, KRAS G12C, and KRAS G12D. The trained models demonstrated strong predictive performance, achieving high correlation coefficients on independent test sets. To further validate the predictive capability of the models, two compounds identified as high-confidence candidates, Cobimetinib and Etrasimod, were selected for experimental evaluation. In vitro testing revealed measurable IC50 and E-max values, with both compounds exhibiting preferential activity in KRAS G12D cellular backgrounds. While additional biochemical and pathway-level studies are required to confirm direct target engagement, these results support the model's utility in prioritizing candidate compounds with allele-specific activity profiles. Overall, this study provides a data-driven framework for identifying potential KRAS-targeted therapies and highlights the value of integrating machine learning predictions with experimental validation.

The work of Martin Karplus, who passed away on Dec 28 2024, was at the forefront of computational chemistry and molecular biophysics for a period of more than sixty years. His career started with a PhD in theoretical chemistry at Caltech with Linus Pauling in 1953 . After performing leading research in molecular quantum chemistry for two decades, in the 1970s he began to incorporate work on biological systems, and his work was instrumental in creating modern computational molecular biophysics. In 2013, he was awarded the Nobel Prize together with Arieh Warshel and Michael Levitt “for the development of multiscale models for complex chemical systems”. This article aims at briefly reviewing the main achievements of the research performed in the Karplus lab from the point of view of the people working under his unique mentorship.

Phosphates are essential for life in all organisms, playing key roles in nucleic acids, signaling, energy transfer, and biosynthesis. We conducted a quantitative analysis of nucleoside triphosphate (NTP) processing enzymes across all enzymatic reactions, revealing their dominance in phosphate reactivity with ATP as the most prevalent substrate. Two main reaction types occur predominantly: cleavage resulting in pyrophosphate release or phosphate release/transfer. The large majority of NTP processing enzymes require divalent Mg2+ ions in a mechanistically analogous manner. However, the catalytically competent metal-ion coordination remains unclear for many NTP-processing enzymes. By examining a vast data set of crystallographic structures, we identified a universal "Mg-pinch" motif, confirming our structural hypothesis for almost all NTP processing enzymes. We hypothesized and subsequently confirmed that the catalytic Mg2+ typically coordinates the two phosphate groups between which the P-O bond is cleaved. We present a comprehensive analysis of NTP processing superfamilies across all species, determining distinct enzyme active site structures. We highlight exceptional cases and propose challenging superfamilies that lack sufficient structural data to determine precise active site coordination. DFT-based QM/MM calculations with full electrostatic embedding support a mechanistic interpretation in which appropriately coordinated Mg2+ ions contribute to catalysis by electrostatic preorganization and polarization of the reacting phosphate groups. The Mg-pinch motif provides a mechanistic framework for understanding the catalytic role of metal ions in NTP processing. Our findings offer insights into enzyme evolution, provide a basis for rational enzyme engineering, and could inform the development of novel therapeutics targeting NTP processing enzymes.

Computing electrostatic interactions remains the bottleneck of molecular dynamics (MD) simulations despite more than a century of effort in developing methods to accelerate the calculation. Previously, we have developed the spherical grids and treecode and Gauss-Legendre-spherical-t (GLST) algorithms for electrostatic interactions. Here, we explain the computational details and discuss the performance of GLST. The GLST algorithm achieves O(N) scaling and should be less demanding in parallel communication compared with the widely used particle mesh Ewald method and likely comparable to the communication costs of the fast multipole method. We find that GLST is suitable for rapid calculation of long-range electrostatic interactions in MD simulations as it has highly tunable accuracy and should scale well on massively parallel computing architectures. The GLST software presented here is available as a standalone library on GitHub.

Iron-sulfur (Fe-S) clusters are critical cofactors in metalloproteins, essential for cellular processes such as energy production, DNA repair, enzymatic catalysis, and metabolic regulation. While Fe-S cluster functions are intimately linked to their redox properties, their precise roles in many proteins remain unclear. In this study, we present a regression model based on experimental redox potential (E m ) data, utilizing only two features: the Fe-S cluster's total charge and the Fe atoms' average valence. This model achieves a high correlation with experimental data (R 2 = 0.82) and an average prediction error of 0.12 V. Applying this model across the Protein Data Bank, we predict E m values for all cataloged Fe-S clusters, uncovering redox potential trends across diverse cluster classes. The computed redox potentials showed strong agreement with experimental values, achieving an overall accuracy of 88%. This streamlined, computationally accessible approach enhances the annotation and mechanistic understanding of Fe-S proteins, offering new insights into the redox variability of electron transport proteins. Our model holds promise for advancing studies of metalloprotein function and facilitating the design of bioinspired redox systems.

The thermodynamics of arginine-phosphate binding is key to cellular signaling, protein-nucleic acid interactions, and membrane protein dynamics. In biomolecules, monoester phosphates are typically employed as strong electrostatic anchors strategically placed in switch domains to mediate specific interactions. In the diester configuration, phosphate groups act as ubiquitous connectors in all nucleic acids and polar lipids, while also engaging in less specific but multiple electrostatic interactions. Here, we employ isothermal titration calorimetry and a set of small-molecule models and peptides to benchmark the ability of the CHARMM force field to accurately reproduce these interactions. We observe good agreement between isothermal titration calorimetry and computational results for methylguanidinium (MGUA) with glycerol and glucose phosphates (MGUA-Gly3P, MGUA-Glu6P), and for arginine-glycine-arginine peptide with inositol triphosphate (RGR-IP3) systems, with experimental binding energies of -3.30 ± 0.30, -3.89 ± 0.30, and -8.96 ± 0.17 kcal/mol, compared with computational values of -4.08 ± 0.00, -4.20 ± 0.00, and -9.17 ± 0.20 kcal/mol, respectively. However, the experimental binding energy of -2.24 ± 0.71 kcal/mol between MGUA and dimethylphosphate in a diester configuration was significantly underestimated in CHARMM computations (-0.51 ± 0.01 kcal/mol). The force field was, therefore, refined by reducing the Lennard-Jones Rmin parameter from 3.55 to 3.405 Å for a specific interaction involving nitrogen and oxygen atoms in MGUA-dimethylphosphate. Our study brings another experimental means for fine-tuning force field parameters for the phosphates in two distinct configurations and enhances the accuracy of modeling nucleic acids, lipids, and membrane proteins.

Glycosaminoglycans (GAGs) are long, anionic polysaccharides abundant in the extracellular matrix and lysosomes, where their electrostatic interactions with proteins are essential for biological function. Computational studies of GAG-containing systems remain challenging due to their significant charge density and conformational flexibility. Here we benchmark two widely used force fields, ff14SB/GLYCAM06 and CHARMM36m, for three experimentally characterized protein–GAG complexes. Both approaches reproduce the general structural and energetic features of GAG–protein interactions. ff14SB/GLYCAM06 yields highly stable trajectories and consistent energetic profiles, whereas CHARMM36m more accurately captures GAG-induced conformational dynamics. These results establish practical benchmarks for accurate atomistic simulations of GAG–protein assemblies and inform future developments in biomolecular force fields.

Here, we use the frequency of the atomic hybridizations (s, sp, sp2, and sp3) of each atom type (H, C, N, O, S, etc.) within a molecule to predict the IC50s of drug-like molecules, focusing on compounds targeting the Thrombin, Estrogen Receptor alpha, and Phosphodiesterase 5A proteins. The Neural Network and Random Forest models yield high correlation coefficients (R2) and low mean square error (MSE) using only 19 features. The atomic hybridizations have been used previously to calculate the molecular polarizability using a simple empirical model (Miller et al. JACS 1979). We show that the atomic hybridizations may also be used to accurately predict the molecular polarizabilities of these molecules. The results show the importance of the induced polarization in protein-ligand binding. Furthermore, the variation in R2 and MSE for the different target proteins indicates that the contribution of the induced polarization to the binding energies is different for different target proteins.

We present four tree-based machine learning models for protein pKa prediction. The four models, Random Forest, Extra Trees, eXtreme Gradient Boosting (XGBoost) and Light Gradient Boosting Machine (LightGBM), were trained on three experimental PDB and pKa datasets, two of which included a notable portion of internal residues. We observed similar performance among the four machine learning algorithms. The best model trained on the largest dataset performs 37% better than the widely used empirical pKa prediction tool PROPKA. The overall RMSE for this model is 0.69, with surface and buried RMSE values being 0.56 and 0.78, respectively, considering six residue types (Asp, Glu, His, Lys, Cys and Tyr), and 0.63 when considering Asp, Glu, His and Lys only. We provide pKa predictions for proteins in human proteome from the AlphaFold Protein Structure Database and observed that 1% of Asp/Glu/Lys residues have highly shifted pKa values close to the physiological pH.

Induced polarization plays a pivotal role in ligand-protein binding by enhancing both the specificity and strength of molecular interactions. As a ligand approaches a protein, their respective electronic clouds redistribute in response to their electrostatic fields a phenomenon governed by induced polarization. The response of a molecules electron density to an external field is quantitatively described by its polarizability tensor. In this study, we calculated polarizability tensors for thousands of drug-like molecules from the CHEMBL database, focusing on compounds targeting the Thrombin, Estrogen Receptor alpha, and Phosphodiesterase 5A proteins using Density Functional Theory (DFT). We show that a machine learning model based on atomic hybridization accurately predict the polarizabilities eigenvalues, calculated with DFT. Then, we build a neural network and random forest models to predict the IC50s based on the same features. The success of these models, despite utilizing a limited number of features, underscores the critical role of induced polarizabilities in determining binding energies. ### Competing Interest Statement The authors have declared no competing interest.

We propose a novel classification algorithm, the Boltzmann Classifier, inspired by the thermodynamic principles underlying the Boltzmann distribution. Our method computes a probabilistic estimate for each class based on an energy function derived from feature-wise deviations between input samples and class-specific centroids. The resulting probabilities are proportional to the exponential negative energies, normalized across classes, analogous to the Boltzmann distribution used in statistical mechanics. In addition, the KT variable can be used to allow the high energy states to be more accessible, which allows the tuning of their probabilities as needed. We evaluate the model performance on several datasets from different applications. The model achieves a high accuracy, which indicates that the Boltzmann Classifier is competitive with standard models like logistic regression and k-nearest neighbors while offering a thermodynamically motivated probabilistic interpretation. our classifier does not require iterative optimization or backpropagation and is thus computationally efficient and easy to integrate into existing workflows. This work demonstrates how ideas from physics can inform new directions in machine learning, providing a foundation for interpretable, energy-based decision-making systems.