Polarizable Force Fields for Biomolecular Simulations: Recent Advances and Applications

Absolute and Relative Binding Free Energy Calculations of Nucleotides to Multiple Protein Classes

Polyphosphate nucleotides, such as ATP, ADP, GTP, and GDP, play a crucial role in modulating protein functions through binding and/or catalytically activating proteins (enzymes). However, accurately calculating the binding free energies for these charged and flexible ligands poses challenges due to slow conformational relaxation and the limitations of force fields. In this study, we examine the accuracy and reliability of alchemical free energy simulations with fixed-charge force fields for the binding of four nucleotides to nine proteins of various classes, including kinases, ATPases, and GTPases. Our results indicate that the alchemical simulations effectively reproduce experimental binding free energies for all proteins that do not undergo significant conformational changes between their triphosphate nucleotide-bound and diphosphate nucleotide-bound states, with 87.5% (7 out of 8) of the absolute binding free energy results for 4 proteins within ±2 kcal/mol of experimental values and 88.9% (8 out of 9) of the relative binding free energy results for 9 proteins within ±3 kcal/mol of experimental values. However, our calculations show significant inaccuracies when divalent ions are included, suggesting that nonpolarizable force fields may not accurately capture interactions involving these ions. Additionally, the presence of highly charged and flexible ligands necessitates extensive conformational sampling to account for the long relaxation times associated with long-range electrostatic interactions. The simulation strategy presented here, along with its demonstrated accuracy across multiple protein classes, will be valuable for predicting the binding of nucleotides or their analogs to protein targets.

Do selectivity filter carbonyls in K+ channels flip away from the pore? Two-dimensional infrared spectroscopy study

Molecular dynamics simulations revealed that the carbonyls of the Val residue in the conserved selectivity filter sequence TVGTG of potassium ion channels can flip away from the pore to form hydrogen bonds with the network of water molecules residing behind the selectivity filter. Such a configuration has been proposed to be relevant for C-type inactivation. Experimentally, X-ray crystallography of the KcsA channel admits the possibility that the Val carbonyls can flip, but it cannot decisively confirm the existence of such a configuration. In this study, we combined molecular dynamics simulations and line shape theory to design two-dimensional infrared spectroscopy experiments that can corroborate the existence of the selectivity filter configuration with flipped Val carbonyls. This ability to distinguish between flipped and unflipped carbonyls is based on the varying strength of the electric field inside and outside the pore, which is directly linked to carbonyl stretching frequencies that can be resolved using infrared spectroscopy.

Theoretically grounded approaches to account for polarization effects in fixed-charge force fields

Non-polarizable, or fixed-charge, force fields are the workhorses of most molecular simulation studies. They attempt to describe the potential energy surface (PES) of the system by including polarization effects in an implicit way. This has historically been done in a rather empirical and ad hoc manner. Recent theoretical treatments of polarization, however, offer promise for getting the most out of fixed-charge force fields by judicious choice of parameters (most significantly the net charge or dipole moment of the model) and application of post facto polarization corrections. This Perspective describes these polarization theories, namely the "halfway-charge" theory and the molecular dynamics in electronic continuum theory, and shows that they lead to qualitatively (and often, quantitatively) similar predictions. Moreover, they can be reconciled into a unified approach to construct a force field development workflow that can yield non-polarizable models with charge/dipole values that provide an optimal description of the PES. Several applications of this approach are reviewed, and avenues for future research are proposed.

Force Fields for Deep Eutectic Mixtures: Application to Structure, Thermodynamics and 2D-Infrared Spectroscopy

Parametrizing energy functions for ionic systems can be challenging. Here, the total energy function for an eutectic system consisting of water, SCN-, K+ and acetamide is improved vis-a-vis experimentally measured properties. Given the importance of electrostatic interactions, two different types of models are considered: the first (model M0) uses atom-centered multipole whereas the other two (models M1 and M2) are based on fluctuating minimal distributed charges (fMDCM) that respond to geometrical changes of SCN-. The Lennard-Jones parameters of the anion are adjusted to best reproduce experimentally known hydration free energies and densities, which are matched to within a few percent for the final models irrespective of the electrostatic model. Molecular dynamics simulations of the eutectic mixtures with varying water content (between 0 and 100%) yield radial distribution functions and frequency correlation functions for the CN-stretch vibration. Comparison with experiments indicates that models based on fMDCM are considerably more consistent than those using multipoles. Computed viscosities from models M1 and M2 are within 30% of measured values and their change with increasing water content is consistent with experiments. This is not the case for model M0.

Exploring Ion Polarizabilities and Their Correlation with van der Waals Radii: A Theoretical Investigation

Polarizability (α) is a fundamental property which measures the tendency of the electron cloud of an atom, ion, or molecule to be distorted by electric field. Polarizability contributes to important physical properties such as molecular interactions or dielectric constants; thus, it is essential to have accurate polarizabilities in molecular simulations. However, it remains a challenge to develop polarizable force fields (FFs) for ions in computational chemistry. In particular, a comprehensive set of polarizabilities for ions has not been derived. Herein, we derived a systematic set of polarizabilities for atoms and ions across the periodic table based on high-level quantum mechanics calculations. These values have excellent agreement with experimental data. Furthermore, we examined the relationship between the obtained polarizabilities and the van der Waals (VDW) radii (RVDW) that we previously determined (J. Chem. Theory Comput., 2023, 19, 2064). Two relationships, RVDW ∝ α1/7 and RVDW ∝ α1/3, proposed in previous studies were examined in the present work. Our results indicated the former relationship, which was derived based on the quantum harmonic oscillator model, prevails for atoms and cations, but neither relationship provides a satisfactory fit for anions. This is consistent with the tight-binding nature of the electrons in atoms and cations, while it is more challenging to quantify the polarizabilities of anions because of their more dispersed electron clouds. Moreover, we compared different approaches to determine the dispersion coefficients, including the London equation, Slater-Kirkwood equation, symmetry-adapted perturbation theory (SAPT) calculations, and time-dependent density functional theory method, along with the approach based on VDW constants. Our results indicated that although different approaches predict deviated magnitudes for the dispersion coefficients, their predictions are highly correlated, implying that each of these approaches can be used to evaluate dispersion interactions after proper scaling. Finally, we have developed a parametrization strategy for the 12-6-4 model based on the obtained insights. We specifically compared the performance of the 12-6-4 model with SAPT and SobEDA analyses to model interactions involving Na+/Mg2+ and various ligands containing He, Ne, Ar, H2O, NH3, [H2PO4]-, and [HPO4]2-. Our results demonstrate that the 12-6-4 parameters effectively reproduce both the total interaction energy and the individual energy components (electrostatics, exchange-repulsion, dispersion, and induction), highlighting the physical robustness of the 12-6-4 model and the effectiveness of our parametrization approach. This study has significant implications for advancing the development of next-generation ion models and polarizable FFs.

Effect of Protein-Polarized Ligand Charges on Relative Protein Ligand Binding Affinities

A major challenge in computer-aided drug design is predicting relative binding energies of different molecules to a target protein using fast and accurate free-energy calculation methods. Free-energy calculations are primarily computed by utilizing classical molecular dynamics simulations based on all-atom force fields (FF) to model the interactions in the system. The present standard classical all-atom FFs contain fixed partial charges on the atoms, and hence electrostatic interactions are modeled between them. The parametrization process to determine these partial charges usually relies on quantum mechanics or semiempirical calculations of the molecule in the gas phase or homogeneous water surrounding. These present standard parametrization schemes of the partial charges neglect, therefore, polarization effects from the protein surrounding. The absence of protein polarization effects can lead to significant errors in free-energy calculations in proteins. We present a parametrization scheme for the partial charges of ligands, named protein-induced polarization (PIP) charges, which account for the electrostatic polarization due to the protein surrounding. The scheme involves single-point quantum mechanics/molecular mechanics calculations of the ligand charges in the protein/water surrounding. Using PIP ligand partial charges, we have calculated the relative binding free energies (RBFEs) of well-studied protein-ligand systems. We show here that RBFEs computed with PIP charges are either significantly improved or at least comparable to those computed with nonpolarized standard GAFF charges. Overall, we present a simple-to-use parametrization scheme to include protein polarization in any type of binding free-energy calculations. The parametrization scheme increases the accuracy in RBFE calculations, while it does not add significant computation time to standard parametrization procedures.

Ion-Specific Effects on Ion and Polyelectrolyte Solvation

Ion-specific effects on aqueous solvation of monovalent counter ions, Na+ ${^+ }$ , K+ ${^+ }$ , Cl- ${^- }$ , and Br- ${^- }$ , and two model polyelectrolytes (PEs), poly(styrene sulfonate) (PSS) and poly(diallyldimethylammonium) (PDADMA) were here studied with ab initio molecular dynamics (AIMD) and classical molecular dynamics (MD) simulations based on the OPLS-aa force-field which is an empirical fixed point-charge force-field. Ion-specific binding to the PE charge groups was also characterized. Both computational methods predict similar response for the solvation of the PEs but differ notably in description of ion solvation. Notably, AIMD captures the experimentally observed differences in Cl- ${^- }$ and Br- ${^- }$ anion solvation and binding with the PEs, while the classical MD simulations fail to differentiate the ion species response. Furthermore, the findings show that combining AIMD with the computationally less costly classical MD simulations allows benefiting from both the increased accuracy and statistics reach.

Importance of Polarizable Embedding for Computing Optical Rotation: The Case of Camphor in Ethanol

Solvation effects on optical rotation are notoriously challenging to model for computational chemistry, as the specific rotatory power of a molecule can vary wildly going from apolar to polar or even protic solvents. To address such a problem, we present a polarizable embedding implementation of an electric and magnetic response property based on density functional theory and the AMOEBA polarizable force field, and apply such an implementation to the study of the optical rotation of camphor in ethanol. By comparing a continuum model, and electrostatic and polarizable embedding QM/MM models, we observe that accounting for the environment's polarization gives rise to not only a different quantitative prediction, in very good agreement with experiments for the QM/AMOEBA model, but also to a very different qualitative picture, with the values of the optical rotation computed along a classical molecular dynamics trajectory with electrostatic embedding being statistically uncorrelated to the ones obtained with the polarizable description.

FeNNol: An efficient and flexible library for building force-field-enhanced neural network potentials

Neural network interatomic potentials (NNPs) have recently proven to be powerful tools to accurately model complex molecular systems while bypassing the high numerical cost of ab initio molecular dynamics simulations. In recent years, numerous advances in model architectures as well as the development of hybrid models combining machine-learning (ML) with more traditional, physically motivated, force-field interactions have considerably increased the design space of ML potentials. In this paper, we present FeNNol, a new library for building, training, and running force-field-enhanced neural network potentials. It provides a flexible and modular system for building hybrid models, allowing us to easily combine state-of-the-art embeddings with ML-parameterized physical interaction terms without the need for explicit programming. Furthermore, FeNNol leverages the automatic differentiation and just-in-time compilation features of the Jax Python library to enable fast evaluation of NNPs, shrinking the performance gap between ML potentials and standard force-fields. This is demonstrated with the popular ANI-2x model reaching simulation speeds nearly on par with the AMOEBA polarizable force-field on commodity GPUs (graphics processing units). We hope that FeNNol will facilitate the development and application of new hybrid NNP architectures for a wide range of molecular simulation problems.

Condensed Matter Systems Exposed to Radiation: Multiscale Theory, Simulations, and Experiment

This roadmap reviews the new, highly interdisciplinary research field studying the behavior of condensed matter systems exposed to radiation. The Review highlights several recent advances in the field and provides a roadmap for the development of the field over the next decade. Condensed matter systems exposed to radiation can be inorganic, organic, or biological, finite or infinite, composed of different molecular species or materials, exist in different phases, and operate under different thermodynamic conditions. Many of the key phenomena related to the behavior of irradiated systems are very similar and can be understood based on the same fundamental theoretical principles and computational approaches. The multiscale nature of such phenomena requires the quantitative description of the radiation-induced effects occurring at different spatial and temporal scales, ranging from the atomic to the macroscopic, and the interlinks between such descriptions. The multiscale nature of the effects and the similarity of their manifestation in systems of different origins necessarily bring together different disciplines, such as physics, chemistry, biology, materials science, nanoscience, and biomedical research, demonstrating the numerous interlinks and commonalities between them. This research field is highly relevant to many novel and emerging technologies and medical applications.

Advancing Force Fields Parameterization: A Directed Graph Attention Networks Approach

Force fields (FFs) are an established tool for simulating large and complex molecular systems. However, parametrizing FFs is a challenging and time-consuming task that relies on empirical heuristics, experimental data, and computational data. Recent efforts aim to automate the assignment of FF parameters using pre-existing databases and on-the-fly ab initio data. In this study, we propose a graph-based force field (GB-FFs) model to directly derive parameters for the Generalized Amber Force Field (GAFF) from chemical environments and research into the influence of functional forms. Our end-to-end parametrization approach predicts parameters by aggregating the basic information in directed molecular graphs, eliminating the need for expert-defined procedures and enhances the accuracy and transferability of GAFF across a broader range of molecular complexes. Simulation results are compared to the original GAFF parametrization. In practice, our results demonstrate an improved transferability of the model, showcasing its improved accuracy in modeling intermolecular and torsional interactions, as well as improved solvation free energies. The optimization approach developed in this work is fully applicable to other nonpolarizable FFs as well as to polarizable ones.

In Silico Chemical Experiments in the Age of AI: From Quantum Chemistry to Machine Learning and Back

Computational chemistry is an indispensable tool for understanding molecules and predicting chemical properties. However, traditional computational methods face significant challenges due to the difficulty of solving the Schrödinger equations and the increasing computational cost with the size of the molecular system. In response, there has been a surge of interest in leveraging artificial intelligence (AI) and machine learning (ML) techniques to in silico experiments. Integrating AI and ML into computational chemistry increases the scalability and speed of the exploration of chemical space. However, challenges remain, particularly regarding the reproducibility and transferability of ML models. This review highlights the evolution of ML in learning from, complementing, or replacing traditional computational chemistry for energy and property predictions. Starting from models trained entirely on numerical data, a journey set forth toward the ideal model incorporating or learning the physical laws of quantum mechanics. This paper also reviews existing computational methods and ML models and their intertwining, outlines a roadmap for future research, and identifies areas for improvement and innovation. Ultimately, the goal is to develop AI architectures capable of predicting accurate and transferable solutions to the Schrödinger equation, thereby revolutionizing in silico experiments within chemistry and materials science.

Convergence criteria for single-step free-energy calculations: the relation between the Π bias measure and the sample variance

Free energy calculations play a crucial role in simulating chemical processes, enzymatic reactions, and drug design. However, assessing the reliability and convergence of these calculations remains a challenge. This study focuses on single-step free-energy calculations using thermodynamic perturbation. It explores how the sample distributions influence the estimated results and evaluates the reliability of various convergence criteria, including Kofke's bias measure Π and the standard deviation of the energy difference ΔU, σ ΔU . The findings reveal that for Gaussian distributions, there is a straightforward relationship between Π and σ ΔU , free energies can be accurately approximated using a second-order cumulant expansion, and reliable results are attainable for σ ΔU up to 25 kcal mol-1. However, interpreting non-Gaussian distributions is more complex. If the distribution is skewed towards more positive values than a Gaussian, converging the free energy becomes easier, rendering standard convergence criteria overly stringent. Conversely, distributions that are skewed towards more negative values than a Gaussian present greater challenges in achieving convergence, making standard criteria unreliable. We propose a practical approach to assess the convergence of estimated free energies.

Aromatic and arginine content drives multiphasic condensation of protein-RNA mixtures

Multiphasic architectures are found ubiquitously in biomolecular condensates and are thought to have important implications for the organization of multiple chemical reactions within the same compartment. Many of these multiphasic condensates contain RNA in addition to proteins. Here, we investigate the importance of different interactions in multiphasic condensates comprising two different proteins and RNA using computer simulations with a residue-resolution coarse-grained model of proteins and RNA. We find that in multilayered condensates containing RNA in both phases, protein-RNA interactions dominate, with aromatic residues and arginine forming the key stabilizing interactions. The total aromatic and arginine content of the two proteins must be appreciably different for distinct phases to form, and we show that this difference increases as the system is driven toward greater multiphasicity. Using the trends observed in the different interaction energies of this system, we demonstrate that we can also construct multilayered condensates with RNA preferentially concentrated in one phase. The "rules" identified can thus enable the design of synthetic multiphasic condensates to facilitate further study of their organization and function.

The multi-scale polarizable pseudo-particle solvent coarse-grained approach: From NaCl salt solutions to polyelectrolyte hydration

We discuss key parameters that affect the reliability of hybrid simulations in the aqueous phase based on an efficient multi-scale coarse-grained polarizable pseudo-particle approach, denoted as pppl, to model the solvent water, whereas solutes are modeled using an all atom polarizable force field. Among those parameters, the extension of the solvent domain (SD) at the solute vicinity (domain in which each solvent particle corresponds to a single water molecule) and the magnitude of solute/solvent short range polarization damping effects are shown to be pivotal to model NaCl salty aqueous solutions and the hydration of charged systems, such as the hydrophobic polyelectrolyte polymer that we have recently investigated [Masella et al., J. Chem. Phys. 155, 114903 (2021)]. Strong short range damping is pivotal to simulate aqueous salt NaCl solutions at moderate concentration (up to 1.0M). The SD extension (as well as short range damping) has a weak effect on the polymer conformation; however, it plays a pivotal role in computing accurate polymer/solvent interaction energies. As the pppl approach is up to two orders of magnitude computationally more efficient than all atom polarizable force field methods, our results show it to be an efficient alternative route to investigate the equilibrium properties of complex charged molecular systems in extended chemical environments.

Water Solvent Reorganization upon Ultrafast Resonant Stimulated X-ray Raman Excitation of a Metalloporphyrin Dimer

We propose an X-ray Raman pump-X-ray diffraction probe scheme to follow solvation dynamics upon charge migration in a solute molecule. The X-ray Raman pump selectively prepares a valence electronic wavepacket in the solute, while the probe provides information about the entire molecular ensemble. A combination of molecular dynamics and ab initio quantum chemistry simulations is applied to a Zn-Ni porphyrin dimer in water. Using time-resolved X-ray diffraction and pair distribution functions, we extracted solvation shell dynamics.

Evaluating Polarizable Biomembrane Simulations against Experiments

Owing to the increase of available computational capabilities and the potential for providing a more accurate description, polarizable molecular dynamics force fields are gaining popularity in modeling biomolecular systems. It is, however, crucial to evaluate how much precision is truly gained with increasing cost and complexity of the simulation. Here, we leverage the NMRlipids open collaboration and Databank to assess the performance of available polarizable lipid models─the CHARMM-Drude and the AMOEBA-based parameters─against high-fidelity experimental data and compare them to the top-performing nonpolarizable models. While some improvement in the description of ion binding to membranes is observed in the most recent CHARMM-Drude parameters, and the conformational dynamics of AMOEBA-based parameters are excellent, the best nonpolarizable models tend to outperform their polarizable counterparts for each property we explored. The identified shortcomings range from inaccuracies in describing the conformational space of lipids to excessively slow conformational dynamics. Our results provide valuable insights for the further refinement of polarizable lipid force fields and for selecting the best simulation parameters for specific applications.

A Vision for the Future of Multiscale Modeling

The rise of modern computer science enabled physical chemistry to make enormous progresses in understanding and harnessing natural and artificial phenomena. Nevertheless, despite the advances achieved over past decades, computational resources are still insufficient to thoroughly simulate extended systems from first principles. Indeed, countless biological, catalytic and photophysical processes require ab initio treatments to be properly described, but the breadth of length and time scales involved makes it practically unfeasible. A way to address these issues is to couple theories and algorithms working at different scales by dividing the system into domains treated at different levels of approximation, ranging from quantum mechanics to classical molecular dynamics, even including continuum electrodynamics. This approach is known as multiscale modeling and its use over the past 60 years has led to remarkable results. Considering the rapid advances in theory, algorithm design, and computing power, we believe multiscale modeling will massively grow into a dominant research methodology in the forthcoming years. Hereby we describe the main approaches developed within its realm, highlighting their achievements and current drawbacks, eventually proposing a plausible direction for future developments considering also the emergence of new computational techniques such as machine learning and quantum computing. We then discuss how advanced multiscale modeling methods could be exploited to address critical scientific challenges, focusing on the simulation of complex light-harvesting processes, such as natural photosynthesis. While doing so, we suggest a cutting-edge computational paradigm consisting in performing simultaneous multiscale calculations on a system allowing the various domains, treated with appropriate accuracy, to move and extend while they properly interact with each other. Although this vision is very ambitious, we believe the quick development of computer science will lead to both massive improvements and widespread use of these techniques, resulting in enormous progresses in physical chemistry and, eventually, in our society.

Anti-Aging Potential of the Two Major Flavonoids Occurring in Asian Water Lily Using In Vitro and In Silico Molecular Modeling Assessments

Our previous study investigated the major flavonoids and antioxidant potential of Asian water lily (Nymphaea lotus L., family Nymphaeaceae) stamens and perianth extracts. Quercetin-3-O-rhamnoside (Que-3-Rha) and kaempferol-3-O-galactoside (Kae-3-Gal) were reported as the two most prominent flavonoids found in these extracts. Many flavonoids have been reported on the skin anti-aging effect that are useful for cosmeceutical/phytopharmaceutical application. However, Que-3-Rha and Kae-3-Gal occurring in this medicinal plant have not yet been evaluated for their ability to inhibit skin-aging enzymes. Therefore, this study aimed (1) to assess the enzyme inhibitory activity of Que-3-Rha and Kae-3-Gal, and (2) to conduct molecular modeling of these compounds against critical enzymes involved in skin aging such as collagenase, elastase, and tyrosinase. In vitro enzymatic assays demonstrated that both of the two most prominent flavonoids exhibited moderate to good inhibitory activity toward these enzymes. These experimental findings were supported by molecular docking analysis, which indicated that Que-3-Rha and Kae-3-Gal showed superior binding affinity to the target enzymes compared to the positive controls. Additionally, computational predictions suggested favorable skin permeability and no severe toxicity for both compounds. The results from molecular dynamic (MD) simulation revealed that all the complexes remained stable during the 200 ns MD simulation. Structural analyses and binding free energy calculations also supported the inhibitory potential of these two flavonoids against skin-aging enzymes. In conclusion, this study provides valuable insights into the anti-aging potential of the two major flavonoids occurring in this medicinal plant, paving the way for further development of cosmeceutical/phytopharmaceutical products targeting skin aging.

QM/MM-Based Energy Decomposition Analysis Method for Large Systems

In this work, a QM/MM-based EDA method, called GKS-EDA(QM/MM), is proposed. As an extension of GKS-EDA, this scheme divides the total interaction energy into electrostatic, exchange-repulsion, polarization, and correlation/dispersion terms. GKS-EDA(QM/MM) can be applied to describe the interactions of large-scale systems combined with various QM/MM platforms. By using the examples of a hydrated hydronium ion complex in water solution, the barnase-barstar complex, and MMP-13-pyrimidinetrione in a metalloprotein, the capability of GKS-EDA(QM/MM) for various interactions in large systems is validated.

Origin of Enantioselectivity in Engineered Cytochrome c-Catalyzed Carbon-Radical FePP Hydrolysis Revealed Using QM/MM (ABEEM Polarizable Force Field) and MD Simulations

The origin of highly efficient asymmetric aminohydroxylation of styrene catalyzed by engineered cytochrome c is investigated by the developed Atom-Bond Electronegativity Equalization Method polarizable force field (ABEEM PFF), which is a combined outcome of electronic and steric effects. Model molecules were used to establish the charge parameters of the ABEEM PFF, for which the bond-stretching and angle-bending parameters were obtained by using a combination of modified Seminario and scan methods. The interactions between carbon-radical Fe-porphyrin (FePP) and waters are simulated by molecular dynamics, which shows a clear preference for the pre-R over the pre-S. This preference is attributed to the hydrogen-bond between the mutated 100S and 101P residues as well as van der Waals interactions, enforcing a specific conformation of the carbon-radical FePP complex within the binding pocket. Meanwhile, the hydrogen-bond between water and the nitrogen atom in the active intermediate dictates the stereochemical outcome. Quantum mechanics/molecular mechanics (QM/MM (ABEEM PFF)) and free-energy perturbation calculations elucidate that the 3RTS is characterized by sandwich-like structure among adjacent amino acid residues, which exhibits greater stability than crowed arrangement in 3STS and enables the R enantiomer to form more favorably. Thus, this study provides mechanistic insight into the catalytic reaction of hemoproteins.

Perspective on Integrative Simulations of Bioenergetic Domains

Bioenergetic processes in cells, such as photosynthesis or respiration, integrate many time and length scales, which makes the simulation of energy conversion with a mere single level of theory impossible. Just like the myriad of experimental techniques required to examine each level of organization, an array of overlapping computational techniques is necessary to model energy conversion. Here, a perspective is presented on recent efforts for modeling bioenergetic phenomena with a focus on molecular dynamics simulations and its variants as a primary method. An overview of the various classical, quantum mechanical, enhanced sampling, coarse-grained, Brownian dynamics, and Monte Carlo methods is presented. Example applications discussed include multiscale simulations of membrane-wide electron transport, rate kinetics of ATP turnover from electrochemical gradients, and finally, integrative modeling of the chromatophore, a photosynthetic pseudo-organelle.

Solvation of Model Biomolecules in Choline-Aminoate Ionic Liquids: A Computational Simulation Using Polarizable Force Fields

One can foresee a very near future where ionic liquids will be used in applications such as biomolecular chemistry or medicine. The molecular details of their interaction with biological matter, however, are difficult to investigate due to the vast number of combinations of both the biological systems and the variety of possible liquids. Here, we provide a computational study aimed at understanding the interaction of a special class of biocompatible ionic liquids (choline-aminoate) with two model biological systems: an oligopeptide and an oligonucleotide. We employed molecular dynamics with a polarizable force field. Our results are in line with previous experimental and computational evidence on analogous systems and show how these biocompatible ionic liquids, in their pure form, act as gentle solvents for protein structures while simultaneously destabilizing DNA structure.

Influence of ion and hydration atmospheres on RNA structure and dynamics: insights from advanced theoretical and computational methods

RNA, a highly charged biopolymer composed of negatively charged phosphate groups, defies electrostatic repulsion to adopt well-defined, compact structures. Hence, the presence of positively charged metal ions is crucial not only for RNA's charge neutralization, but they also coherently decorate the ion atmosphere of RNA to stabilize its compact fold. This feature article elucidates various modes of close RNA-ion interactions, with a special emphasis on Mg2+ as an outer-sphere and inner-sphere ion. Through examples, we highlight how inner-sphere chelated Mg2+ stabilizes RNA pseudoknots, while outer-sphere ions can also exert long-range electrostatic interactions, inducing groove narrowing, coaxial helical stacking, and RNA ring formation. In addition to investigating the RNA's ion environment, we note that the RNA's hydration environment is relatively underexplored. Our study delves into its profound interplay with the structural dynamics of RNA, employing state-of-the-art atomistic simulation techniques. Through examples, we illustrate how specific ions and water molecules are associated with RNA functions, leveraging atomistic simulations to identify preferential ion binding and hydration sites. However, understanding their impact(s) on the RNA structure remains challenging due to the involvement of large length and long time scales associated with RNA's dynamic nature. Nevertheless, our contributions and recent advances in coarse-grained simulation techniques offer insights into large-scale structural changes dynamically linked to the RNA ion atmosphere. In this connection, we also review how different cutting-edge computational simulation methods provide a microscopic lens into the influence of ions and hydration on RNA structure and dynamics, elucidating distinct ion atmospheric components and specific hydration layers and their individual and collective impacts.

Incorporating Neural Networks into the AMOEBA Polarizable Force Field

Neural network potentials (NNPs) offer significant promise to bridge the gap between the accuracy of quantum mechanics and the efficiency of molecular mechanics in molecular simulation. Most NNPs rely on the locality assumption that ensures the model's transferability and scalability and thus lack the treatment of long-range interactions, which are essential for molecular systems in the condensed phase. Here we present an integrated hybrid model, AMOEBA+NN, which combines the AMOEBA potential for the short- and long-range noncovalent atomic interactions and an NNP to capture the remaining local covalent contributions. The AMOEBA+NN model was trained on the conformational energy of the ANI-1x data set and tested on several external data sets ranging from small molecules to tetrapeptides. The hybrid model demonstrated substantial improvements over the baseline models in term of accuracy as the molecule size increased, suggesting its potential as a next-generation approach for chemically accurate molecular simulations.

Efficient Sampling of Cavity Hydration in Proteins with Nonequilibrium Grand Canonical Monte Carlo and Polarizable Force Fields

Prediction of the hydration levels of protein cavities and active sites is important to both mechanistic analysis and ligand design. Due to the unique microscopic environment of these buried water molecules, a polarizable model is expected to be crucial for an accurate treatment of protein internal hydration in simulations. Here we adapt a nonequilibrium candidate Monte Carlo approach for conducting grand canonical Monte Carlo simulations with the Drude polarizable force field. The GPU implementation enables the efficient sampling of internal cavity hydration levels in biomolecular systems. We also develop an enhanced sampling approach referred to as B-walking, which satisfies detailed balance and readily combines with grand canonical integration to efficiently calculate quantitative binding free energies of water to protein cavities. Applications of these developments are illustrated in a solvent box and the polar ligand binding site in trypsin. Our simulation results show that including electronic polarization leads to a modest but clear improvement in the description of water position and occupancy compared to the crystal structure. The B-walking approach enhances the range of water sampling in different chemical potential windows and thus improves the accuracy of water binding free energy calculations.

Molecular driving forces for water adsorption in MOF-808: A comparative analysis with UiO-66

Metal-organic frameworks (MOFs), with their unique porous structures and versatile functionality, have emerged as promising materials for the adsorption, separation, and storage of diverse molecular species. In this study, we investigate water adsorption in MOF-808, a prototypical MOF that shares the same secondary building unit (SBU) as UiO-66, and elucidate how differences in topology and connectivity between the two MOFs influence the adsorption mechanism. To this end, molecular dynamics simulations were performed to calculate several thermodynamic and dynamical properties of water in MOF-808 as a function of relative humidity (RH), from the initial adsorption step to full pore filling. At low RH, the μ3-OH groups of the SBUs form hydrogen bonds with the initial water molecules entering the pores, which triggers the filling of these pores before the μ3-OH groups in other pores become engaged in hydrogen bonding with water molecules. Our analyses indicate that the pores of MOF-808 become filled by water sequentially as the RH increases. A similar mechanism has been reported for water adsorption in UiO-66. Despite this similarity, our study highlights distinct thermodynamic properties and framework characteristics that influence the adsorption process differently in MOF-808 and UiO-66.

Structural Insights and Influence of Terahertz Waves in Midinfrared Region on Kv1.2 Channel Selectivity Filter

Potassium ion channels are the structural basis for excitation transmission, heartbeat, and other biological processes. The selectivity filter is a critical structural component of potassium ion channels, whose structure is crucial to realizing their function. As biomolecules vibrate and rotate at frequencies in the terahertz band, potassium ion channels are sensitive to terahertz waves. Therefore, it is worthwhile to investigate how the terahertz wave influences the selectivity filter of the potassium channels. In this study, we investigate the structure of the selectivity filter of Kv1.2 potassium ion channels using molecular dynamics simulations. The effect of an electric field on the channel has been examined at four different resonant frequencies of the carbonyl group in SF: 36.75 37.06, 37.68, and 38.2 THz. As indicated by the results, 376GLY appears to be the critical residue in the selectivity filter of the Kv1.2 channel. Its dihedral angle torsion is detrimental to the channel structural stability and the transmembrane movement of potassium ions. 36.75 THz is the resonance frequency of the carbonyl group of 376GLY. Among all four frequencies explored, the applied terahertz electric field of this frequency has the most significant impact on the channel structure, negatively impacting the channel stability and reducing the ion permeability by 20.2% compared to the absence of fields. In this study, we simulate that terahertz waves in the mid-infrared frequency region can significantly alter the structure and function of potassium ion channels and that the effects of terahertz waves differ greatly based on frequency.

Multiscale computational simulation of pollutant behavior at water interfaces

The investigation of pollutant behavior at water interfaces is critical to understand pollution in aquatic systems. Computational methods allow us to overcome the limitations of experimental analysis, delivering valuable insights into the chemical mechanisms and structural characteristics of pollutant behavior at interfaces across a range of scales, from microscopic to mesoscopic. Quantum mechanics, all-atom molecular dynamics simulations, coarse-grained molecular dynamics simulations, and dissipative particle dynamics simulations represent diverse molecular interaction calculation methods that can effectively model pollutant behavior at environmental interfaces from atomic to mesoscopic scales. These methods provide a rich variety of information on pollutant interactions with water surfaces. This review synthesizes the advancements in applying typical computational methods to the formation, adsorption, binding, and catalytic conversion of pollutants at water interfaces. By drawing on recent advancements, we critically examine the current challenges and offer our perspective on future directions. This review seeks to advance our understanding of computational techniques for elucidating pollutant behavior at water interfaces, a critical aspect of water research.

A Fast, Convenient, Polarizable Electrostatic Model for Molecular Dynamics

We present an efficient polarizable electrostatic model, utilizing typed, atom-centered polarizabilities and the fast direct approximation, designed for efficient use in molecular dynamics (MD) simulations. The model provides two convenient approaches for assigning partial charges in the context of atomic polarizabilities. One is a generalization of RESP, called RESP-dPol, and the other, AM1-BCC-dPol, is an adaptation of the widely used AM1-BCC method. Both are designed to accurately replicate gas-phase quantum mechanical electrostatic potentials. Benchmarks of this polarizable electrostatic model against gas-phase dipole moments, molecular polarizabilities, bulk liquid densities, and static dielectric constants of organic liquids show good agreement with the reference values. Of note, the model yields markedly more accurate dielectric constants of organic liquids, relative to a matched nonpolarizable force field. MD simulations with this method, which is currently parametrized for molecules containing elements C, N, O, and H, run only about 3.6-fold slower than fixed charge force fields, while simulations with the self-consistent mutual polarization average 4.5-fold slower. Our results suggest that RESP-dPol and AM1-BCC-dPol afford improved accuracy relative to fixed charge force fields and are good starting points for developing general, affordable, and transferable polarizable force fields. The software implementing these approaches has been designed to utilize the force field fitting frameworks developed and maintained by the Open Force Field Initiative, setting the stage for further exploration of this approach to polarizable force field development.

Glass Transition Temperatures and Thermal Conductivities of Polybutadiene Crosslinked with Randomly Distributed Sulfur Chains Using Molecular Dynamic Simulation

The thermal conductivities and glass transition temperatures of polybutadiene crosslinked with randomly distributed sulfur chains having different lengths from mono-sulfur (S1) to octa-sulfur (S8) were investigated. The thermal conductivities of the related models as a function of the heat flux autocorrelation function, applying an equilibrium molecular dynamic (EMD) simulation and the Green-Kubo method, were studied for a wide range of temperatures. The influence of the length of sulfur chains, degree of crosslinking, and molar mass of the crosslinker on the glass transition temperature and final values of thermal conductivities were studied. First, the degree of crosslinking is considered constant for the eight simulation models, from mono-sulfur (S1) to octa-sulfur (S8), while the molar mass of the sulfur is increases. The results show that the thermal conductivities of the crosslinked structure decrease with increasing temperature for each model. Moreover, by increasing the lengths of the sulfur chains and the molar weight of the crosslinker, thermal conductivity increases at a constant temperature. The MD simulation demonstrates that the glass transition temperature and density of the crosslinked structure enhance as the length of the sulfur chains and molar mass of the sulfur increase. Second, the molar weight of sulfur is considered constant in these eight models; therefore, the degree of crosslinking decreases with the increase in the lengths of the sulfur chains. The results show that the thermal conductivities of the crosslinked structure decrease with the increase in the temperature for each model. Moreover, by increasing the lengths of sulfur chains and thus decreasing the degree of crosslinking, the trend in changes in thermal conductivities are almost the same for all of these models, so thermal conductivity is constant for a specific temperature. In addition, the glass transition temperature and density of the crosslinked structure decrease.

Thermodynamics of Metal-Acetate Interactions

Metal ions play crucial roles in protein- and ligand-mediated interactions. They not only act as catalysts to facilitate biological processes but are also important as protein structural elements. Accurately predicting metal ion interactions in computational studies has always been a challenge, and various methods have been suggested to improve these interactions. One such method is the 12-6-4 Lennard-Jones (LJ)-type nonbonded model. Using this model, it has been possible to successfully reproduce the experimental properties of metal ions in aqueous solution. The model includes induced dipole interactions typically ignored in the standard 12-6 LJ nonbonded model. In this we expand the applicability of this model to metal ion-carboxylate interactions. Using 12-6-4 parameters that reproduce the solvation free energies of the metal ions leads to an overestimation of metal ion-acetate interactions, thus, prompting us to fine-tune the model to specifically handle the latter. We also show that the standard 12-6 LJ model significantly falls short in reproducing the experimental binding free energy between acetate and 11 metal ions (Ni(II), Mg(II), Cu(II), Zn(II), Co(II), Cu(I), Fe(II), Mn(II), Cd(II), Ca(II), and Ag(I)). In this study, we describe optimized C4 parameters for the 12-6-4 LJ nonbonded model to be used with three widely employed water models (Transferable Intermolecular Potential with 3 Points (TIP3P), Simple Point Charge Extended (SPC/E), and Optimal Point Charge (OPC) water models). These parameters can accurately match the experimental binding free energy between 11 metal ions and acetate. These parameters can be applied to the study of metalloproteins and transition metal ion channels and transporters, as acetate serves as a representative of the negatively charged amino acid side chains from aspartate and glutamate.

A Critical Assessment on Calculating Vibrational Spectra in Nanostructured Materials

Vibrational spectroscopy is an omnipresent spectroscopic technique to characterize functional nanostructured materials such as zeolites, metal-organic frameworks (MOFs), and metal-halide perovskites (MHPs). The resulting experimental spectra are usually complex, with both low-frequency framework modes and high-frequency functional group vibrations. Therefore, theoretically calculated spectra are often an essential element to elucidate the vibrational fingerprint. In principle, there are two possible approaches to calculate vibrational spectra: (i) a static approach that approximates the potential energy surface (PES) as a set of independent harmonic oscillators and (ii) a dynamic approach that explicitly samples the PES around equilibrium by integrating Newton's equations of motions. The dynamic approach considers anharmonic and temperature effects and provides a more genuine representation of materials at true operating conditions; however, such simulations come at a substantially increased computational cost. This is certainly true when forces and energy evaluations are performed at the quantum mechanical level. Molecular dynamics (MD) techniques have become more established within the field of computational chemistry. Yet, for the prediction of infrared (IR) and Raman spectra of nanostructured materials, their usage has been less explored and remain restricted to some isolated successes. Therefore, it is currently not a priori clear which methodology should be used to accurately predict vibrational spectra for a given system. A comprehensive comparative study between various theoretical methods and experimental spectra for a broad set of nanostructured materials is so far lacking. To fill this gap, we herein present a concise overview on which methodology is suited to accurately predict vibrational spectra for a broad range of nanostructured materials and formulate a series of theoretical guidelines to this purpose. To this end, four different case studies are considered, each treating a particular material aspect, namely breathing in flexible MOFs, characterization of defects in the rigid MOF UiO-66, anharmonic vibrations in the metal-halide perovskite CsPbBr3, and guest adsorption on the pores of the zeolite H-SSZ-13. For all four materials, in their guest- and defect-free state and at sufficiently low temperatures, both the static and dynamic approach yield qualitatively similar spectra in agreement with experimental results. When the temperature is increased, the harmonic approximation starts to fail for CsPbBr3 due to the presence of anharmonic phonon modes. Also, the spectroscopic fingerprints of defects and guest species are insufficiently well predicted by a simple harmonic model. Both phenomena flatten the potential energy surface (PES), which facilitates the transitions between metastable states, necessitating dynamic sampling. On the basis of the four case studies treated in this Review, we can propose the following theoretical guidelines to simulate accurate vibrational spectra of functional solid-state materials: (i) For nanostructured crystalline framework materials at low temperature, insights into the lattice dynamics can be obtained using a static approach relying on a few points on the PES and an independent set of harmonic oscillators. (ii) When the material is evaluated at higher temperatures or when additional complexity enters the system, e.g., strong anharmonicity, defects, or guest species, the harmonic regime breaks down and dynamic sampling is required for a correct prediction of the phonon spectrum. These guidelines and their illustrations for prototype material classes can help experimental and theoretical researchers to enhance the knowledge obtained from a lattice dynamics study.

Computational approaches to investigate fluoride binding, selectivity and transport across the membrane

The use of molecular dynamics (MD) simulations to study biomolecular systems has proven reliable in elucidating atomic-level details of structure and function. In this chapter, MD simulations were used to uncover new insights into two phylogenetically unrelated bacterial fluoride (F-) exporters: the CLCF F-/H+ antiporter and the Fluc F- channel. The CLCF antiporter, a member of the broader CLC family, has previously revealed unique stoichiometry, anion-coordinating residues, and the absence of an internal glutamate crucial for proton import in the CLCs. Through MD simulations enhanced with umbrella sampling, we provide insights into the energetics and mechanism of the CLCF transport process, including its selectivity for F- over HF. In contrast, the Fluc F- channel presents a novel architecture as a dual topology dimer, featuring two pores for F- export and a central non-transported sodium ion. Using computational electrophysiology, we simulate the electrochemical gradient necessary for F- export in Fluc and reveal details about the coordination and hydration of both F- and the central sodium ion. The procedures described here delineate the specifics of these advanced techniques and can also be adapted to investigate other membrane protein systems.

Polarizable AMOEBA Model for Simulating Mg2+·Protein·Nucleotide Complexes

Molecular mechanics (MM) simulations have the potential to provide detailed insights into the mechanisms of enzymes that utilize nucleotides as cofactors. In most cases, the activities of these enzymes also require the binding of divalent cations to catalytic sites. However, modeling divalent cations in MM simulations has been challenging. The inclusion of explicit polarization was considered promising, but despite improvements over nonpolarizable force fields and despite the inclusion of "Nonbonded-fix (NB-fix)" corrections, errors in interaction energies of divalent cations with proteins remain large. Importantly, the application of these models fails to reproduce the experimental structural data on Mg2+·Protein·ATP complexes. Focusing on these complexes, here we provide a systematic assessment of the polarizable AMOEBA model and recommend critical changes that substantially improve its predictive performance. Our key results are as follows. We first show that our recent revision of the AMOEBA protein model (AMOEBABIO18-HFC), which contains high field corrections (HFCs) to induced dipoles, dramatically improves Mg2+-protein interaction energies, reducing the mean absolute error (MAE) from 17 to 10 kcal/mol. This further supports the general applicability of AMOEBABIO18-HFC. The inclusion of many-body NB-fix corrections further reduces MAE to 6 kcal/mol, which amounts to less than 2% error. The errors are estimated with respect to vdW-inclusive density functional theory that we benchmark against CCSD(T) calculations and experiments. We also present a new model of ATP with revised polarization parameters to better capture its high field response, as well as new vdW and dihedral parameters. The ATP model accurately predicts experimental Mg2+-ATP binding free energy in the aqueous phase and provides new insights into how Mg2+ associates with ATP. Finally, we show that molecular dynamics (MD) simulations of Mg2+·Kinase·ATP complexes carried out with these improvements lead to a better agreement in global and local catalytic site structures between MD and X-ray crystallography.

Water Exchange from the Buried Binding Sites of Cytochrome P450 Enzymes 1A2, 2D6, and 3A4 Correlates with Conformational Fluctuations

Human cytochrome P450 enzymes (CYPs) are critical for the metabolism of small-molecule pharmaceuticals (drugs). As such, the prediction of drug metabolism by and drug inhibition of CYP activity is an important component of the drug discovery and design process. Relative to the availability of a wide range of experimental atomic-resolution CYP structures, the development of structure-based CYP activity models has been limited. To better characterize the role of CYP conformational fluctuations in CYP activity, we perform multiple microsecond-scale all-atom explicit-solvent molecular dynamics (MD) simulations on three CYP isoforms, 1A2, 2D6, and 3A4, which together account for the majority of CYP-mediated drug metabolism. The MD simulations employ a variety of positional restraints, ranging from keeping all CYP atoms close to their experimentally determined coordinates to allowing full flexibility. We find that, with full flexibility, large fluctuations in the CYP binding sites correlate with efficient water exchange from these buried binding sites. This is especially true for 1A2, which, when restrained to its crystallographic conformation, is unable to exchange water between the binding site and bulk solvent. These findings imply that, in addition to crystal structures, a representative ensemble of conformational states ought to be included when developing structure-based CYP activity models.

A Multiple-Fidelity Method for Accurate Simulation of MoS2 Properties Using JAX-ReaxFF and Neural Network Potentials

Reactive force field (ReaxFF) is a commonly used force field for modeling chemical reactions at the atomic level. Recently, JAX-ReaxFF, combined with automatic differentiation, has been used to efficiently parametrize ReaxFF. However, its analytical formula may lead to inaccurate predictions. While neural network-based potentials (NNPs) trained on density functional theory-labeled data offer a more accurate method, it requires a large amount of training data to be trained from scratch. To overcome these issues, we present a multiple-fidelity method that combines JAX-ReaxFF and NNP and apply the method on MoS2, a promising two-dimensional semiconductor for flexible electronics. By incorporating implicit prior physical information, ReaxFF can serve as a cost-effective way to generate pretraining data, facilitating more accurate simulations of MoS2. Moreover, in the Mo-S-H system, the pretraining strategy can reduce root-mean-square errors of energy by 20%. This approach can be extended to a wide variety of material systems, accelerating their computational research.

Water inside the Selectivity Filter of a K+ Ion Channel: Structural Heterogeneity, Picosecond Dynamics, and Hydrogen Bonding

Water inside biological ion channels regulates the key properties of these proteins, such as selectivity, ion conductance, and gating. In this article, we measure the picosecond spectral diffusion of amide I vibrations of an isotope-labeled KcsA potassium channel using two-dimensional infrared (2D IR) spectroscopy. By combining waiting time (100-2000 fs) 2D IR measurements of the KcsA channel including 13C18O isotope-labeled Val76 and Gly77 residues with molecular dynamics simulations, we elucidated the site-specific dynamics of water and K+ ions inside the selectivity filter of KcsA. We observe inhomogeneous 2D line shapes with extremely slow spectral diffusion. Our simulations quantitatively reproduce the experiments and show that water is the only component with any appreciable dynamics, whereas K+ ions and the protein are essentially static on a picosecond timescale. By analyzing simulated and experimental vibrational frequencies, we find that water in the selectivity filter can be oriented to form hydrogen bonds with adjacent or nonadjacent carbonyl groups with the reorientation timescales being three times slower and comparable to that of water molecules in liquid, respectively. Water molecules can reside in the cavity sufficiently far from carbonyls and behave essentially like "free" gas-phase-like water with fast reorientation times. Remarkably, no interconversion between these configurations was observed on a picosecond timescale. These dynamics are in stark contrast with liquid water, which remains highly dynamic even in the presence of ions at high concentrations.

Improving the Accuracy of Physics-Based Hydration-Free Energy Predictions by Machine Learning the Remaining Error Relative to the Experiment

The accuracy of computational models of water is key to atomistic simulations of biomolecules. We propose a computationally efficient way to improve the accuracy of the prediction of hydration-free energies (HFEs) of small molecules: the remaining errors of the physics-based models relative to the experiment are predicted and mitigated by machine learning (ML) as a postprocessing step. Specifically, the trained graph convolutional neural network attempts to identify the "blind spots" in the physics-based model predictions, where the complex physics of aqueous solvation is poorly accounted for, and partially corrects for them. The strategy is explored for five classical solvent models representing various accuracy/speed trade-offs, from the fast analytical generalized Born (GB) to the popular TIP3P explicit solvent model; experimental HFEs of small neutral molecules from the FreeSolv set are used for the training and testing. For all of the models, the ML correction reduces the resulting root-mean-square error relative to the experiment for HFEs of small molecules, without significant overfitting and with negligible computational overhead. For example, on the test set, the relative accuracy improvement is 47% for the fast analytical GB, making it, after the ML correction, almost as accurate as uncorrected TIP3P. For the TIP3P model, the accuracy improvement is about 39%, bringing the ML-corrected model's accuracy below the 1 kcal/mol threshold. In general, the relative benefit of the ML corrections is smaller for more accurate physics-based models, reaching the lower limit of about 20% relative accuracy gain compared with that of the physics-based treatment alone. The proposed strategy of using ML to learn the remaining error of physics-based models offers a distinct advantage over training ML alone directly on reference HFEs: it preserves the correct overall trend, even well outside of the training set.

Molecular Docking of Intrinsically Disordered Proteins: Challenges and Strategies

Intrinsically disordered proteins (IDPs) are a novel class of proteins that have established a significant importance and attention within a very short period of time. These proteins are essentially characterized by their inherent structural disorder, encoded mainly by their amino acid sequences. The profound abundance of IDPs and intrinsically disordered regions (IDRs) in the biological world delineates their deep-rooted functionality. IDPs and IDRs convey such extensive functionality through their unique dynamic nature, which enables them to carry out huge number of multifaceted biomolecular interactions and make them "interaction hub" of the cellular systems. Additionally, with such widespread functions, their misfunctioning is also intimately associated with multiple diseases. Thus, understanding the dynamic heterogeneity of various IDPs along with their interactions with respective binding partners is an important field with immense potentials in biomolecular research. In this context, molecular docking-based computational approaches have proven to be remarkable in case of ordered proteins. Molecular docking methods essentially model the biomolecular interactions in both structural and energetic terms and use this information to characterize the putative interactions between the two participant molecules. However, direct applications of the conventional docking methods to study IDPs are largely limited by their structural heterogeneity and demands for unique IDP-centric strategies. Thus, in this chapter, we have presented an overview of current methodologies for successful docking operations involving IDPs and IDRs. These specialized methods majorly include the ensemble-based and fragment-based approaches with their own benefits and limitations. More recently, artificial intelligence and machine learning-assisted approaches are also used to significantly reduce the complexity and computational burden associated with various docking applications. Thus, this chapter aims to provide a comprehensive summary of major challenges and recent advancements of molecular docking approaches in the IDP field for their better utilization and greater applicability.Asp (D).

Force Field Analysis Software and Tools (FFAST): Assessing Machine Learning Force Fields under the Microscope

As the sophistication of machine learning force fields (MLFF) increases to match the complexity of extended molecules and materials, so does the need for tools to properly analyze and assess the practical performance of MLFFs. To go beyond average error metrics and into a complete picture of a model's applicability and limitations, we developed FFAST (force field analysis software and tools): a cross-platform software package designed to gain detailed insights into a model's performance and limitations, complete with an easy-to-use graphical user interface. The software allows the user to gauge the performance of any molecular force field,─such as popular state-of-the-art MLFF models, ─ on various popular data set types, providing general prediction error overviews, outlier detection mechanisms, atom-projected errors, and more. It has a 3D visualizer to find and picture problematic configurations, atoms, or clusters in a large data set. In this paper, the example of the MACE and NequIP models is used on two data sets of interest [stachyose and docosahexaenoic acid (DHA)]─to illustrate the use cases of the software. With this, it was found that carbons and oxygens involved in or near glycosidic bonds inside the stachyose molecule present increased prediction errors. In addition, prediction errors on DHA rise as the molecule folds, especially for the carboxylic group at the edge of the molecule. We emphasize the need for a systematic assessment of MLFF models for ensuring their successful application to the study of dynamics of molecules and materials.

Assessment of Halogen Off-Center Point-Charge Models Using Explicit Solvent Simulations

Compounds containing halogens can form halogen bonds (XBs) with biological targets such as proteins and membranes due to their anisotropic electrostatic potential. To accurately describe this anisotropy, off-center point-charge (EP) models are commonly used in force field methods, allowing the description of XBs at the molecular mechanics and molecular dynamics level. Various EP implementations have been documented in the literature, and despite being efficient in reproducing protein-ligand geometries and sampling of XBs, it is unclear how well these EP models predict experimental properties such as hydration free energies (ΔGhyd), which are often used to validate force field performance. In this work, we report the first assessment of three EP models using alchemical free energy calculations to predict ΔGhyd values. We show that describing the halogen anisotropy using some EP models can lead to a slight improvement in the prediction of the ΔGhyd when compared with the models without EP, especially for the chlorinated compounds; however, this improvement is not related to the establishment of XBs but is most likely due to the improvement of the sampling of hydrogen bonds. We also highlight the importance of the choice of the EP model, especially for the iodinated molecules, since a slight tendency to improve the prediction is observed for compounds with a larger σ-hole but significantly worse results were obtained for compounds that are weaker XB donors.

In-Cell Dynamics: The Next Focus of All-Atom Simulations

The cell is a crowded space where large biomolecules and metabolites are in continuous motion. Great strides have been made in in vitro studies of protein dynamics, folding, and protein-protein interactions, and much new data are emerging of how they differ in the cell. In this Perspective, we highlight the current progress in atomistic modeling of in-cell environments, both bacteria and mammals, with emphasis on classical all-atom molecular dynamics simulations. These simulations have been recently used to capture and characterize functional and non-functional protein-protein interactions, protein folding dynamics of small proteins with varied topologies, and dynamics of metabolites. We further discuss the challenges and efforts for updating modern force fields critical to the progress of cellular environment simulations. We also briefly summarize developments in relevant state-of-the-art experimental techniques. As computational and experimental methodologies continue to progress and produce more directly comparable data, we are poised to capture the complex atomistic picture of the cell.

QMrebind: incorporating quantum mechanical force field reparameterization at the ligand binding site for improved drug-target kinetics through milestoning simulations

Understanding the interaction of ligands with biomolecules is an integral component of drug discovery and development. Challenges for computing thermodynamic and kinetic quantities for pharmaceutically relevant receptor-ligand complexes include the size and flexibility of the ligands, large-scale conformational rearrangements of the receptor, accurate force field parameters, simulation efficiency, and sufficient sampling associated with rare events. Our recently developed multiscale milestoning simulation approach, SEEKR2 (Simulation Enabled Estimation of Kinetic Rates v.2), has demonstrated success in predicting unbinding (koff) kinetics by employing molecular dynamics (MD) simulations in regions closer to the binding site. The MD region is further subdivided into smaller Voronoi tessellations to improve the simulation efficiency and parallelization. To date, all MD simulations are run using general molecular mechanics (MM) force fields. The accuracy of calculations can be further improved by incorporating quantum mechanical (QM) methods into generating system-specific force fields through reparameterizing ligand partial charges in the bound state. The force field reparameterization process modifies the potential energy landscape of the bimolecular complex, enabling a more accurate representation of the intermolecular interactions and polarization effects at the bound state. We present QMrebind (Quantum Mechanical force field reparameterization at the receptor-ligand binding site), an ORCA-based software that facilitates reparameterizing the potential energy function within the phase space representing the bound state in a receptor-ligand complex. With SEEKR2 koff estimates and experimentally determined kinetic rates, we compare and interpret the receptor-ligand unbinding kinetics obtained using the newly reparameterized force fields for model host-guest systems and HSP90-inhibitor complexes. This method provides an opportunity to achieve higher accuracy in predicting receptor-ligand koff rate constants.

Force-field-enhanced neural network interactions: from local equivariant embedding to atom-in-molecule properties and long-range effects

We introduce FENNIX (Force-Field-Enhanced Neural Network InteraXions), a hybrid approach between machine-learning and force-fields. We leverage state-of-the-art equivariant neural networks to predict local energy contributions and multiple atom-in-molecule properties that are then used as geometry-dependent parameters for physically-motivated energy terms which account for long-range electrostatics and dispersion. Using high-accuracy ab initio data (small organic molecules/dimers), we trained a first version of the model. Exhibiting accurate gas-phase energy predictions, FENNIX is transferable to the condensed phase. It is able to produce stable Molecular Dynamics simulations, including nuclear quantum effects, for water predicting accurate liquid properties. The extrapolating power of the hybrid physically-driven machine learning FENNIX approach is exemplified by computing: (i) the solvated alanine dipeptide free energy landscape; (ii) the reactive dissociation of small molecules.

Structure and Dynamics of Confined Water Inside Diphenylalanine Peptide Nanotubes

Diphenylalanine (FF) peptides exhibit a unique ability to self-assemble into nanotubes with confined water molecules playing pivotal roles in their structure and function. This study investigates the structure and dynamics of diphenylalanine peptide nanotubes (FFPNTs) using all-atom molecular dynamics (MD) and grand canonical Monte Carlo combined with MD (GCMC/MD) simulations with both the CHARMM additive and Drude polarizable force fields. The occupancy and dynamics of confined water molecules were also examined. It was found that less than 2 confined water molecules per FF help stabilize the FFPNTs on the x-y plane. Analyses of the kinetics of confined water molecules revealed distinctive transport behaviors for bound and free water, and their respective diffusion coefficients were compared. Our results validate the importance of polarizable force field models in studying peptide nanotubes and provide insights into our understanding of nanoconfined water.

Quantum mechanics based non-bonded force field functions for use in molecular dynamics simulations of materials and systems: The nitrogen and oxygen columns

Accurate Force Fields (FFs) are essential for Molecular Dynamics (MD) simulations of the dynamics of realistic materials in terms of atomic-level interactions. The FF parameters of short-range valence interactions can be derived through Quantum Mechanical (QM) calculations on model systems practical for QM (<300 atoms). Similarly, the dynamic electrostatic interactions can be described with methods such as QEq or PQEq that allow charges and polarization to adjust dynamically. However, accurately extracting long-range van der Waals (vdW) interactions from QM calculations poses challenges due to the absence of a definitive method to distinguish between the different energetic components of electrostatics, polarization, vdW, hydrogen bonding, and valence interactions. To do this we use the Perdew-Burke-Ernzerhof flavor of Density Functional Theory, including empirical D3 vdW corrections, to predict the Equation of State for each element (keeping any covalent bonds fixed), from which we obtain the two-body vdW nonbond potential. Here, we extend these calculations to include non-bonded parameters for the N and O columns of the periodic table so that we now describe columns 15 (N), 16 (O), 17 (F), and 18 (Ne) of the periodic table. For these 20 elements, we find that the two-body vdW potentials can all be mapped to a single universal two-body curve, with just three scaling parameters: Re, De, and L. We refer to this as the Universal NonBond (UNB) potential. We expect this to be useful for new MD simulations and a helpful starting point to obtain UNB parameters for the remainder of the periodic table.

Advances in Computational Approaches for Estimating Passive Permeability in Drug Discovery

Passive permeation of cellular membranes is a key feature of many therapeutics. The relevance of passive permeability spans all biological systems as they all employ biomembranes for compartmentalization. A variety of computational techniques are currently utilized and under active development to facilitate the characterization of passive permeability. These methods include lipophilicity relations, molecular dynamics simulations, and machine learning, which vary in accuracy, complexity, and computational cost. This review briefly introduces the underlying theories, such as the prominent inhomogeneous solubility diffusion model, and covers a number of recent applications. Various machine-learning applications, which have demonstrated good potential for high-volume, data-driven permeability predictions, are also discussed. Due to the confluence of novel computational methods and next-generation exascale computers, we anticipate an exciting future for computationally driven permeability predictions.

LVC/MM: A Hybrid Linear Vibronic Coupling/Molecular Mechanics Model with Distributed Multipole-Based Electrostatic Embedding for Highly Efficient Surface Hopping Dynamics in Solution

We present a theoretical framework for a hybrid linear vibronic coupling model electrostatically embedded into a molecular mechanics environment, termed the linear vibronic coupling/molecular mechanics (LVC/MM) method, for the surface hopping including arbitrary coupling (SHARC) molecular dynamics package. Electrostatic embedding is realized through the computation of interactions between environment point charges and distributed multipole expansions (DMEs, up to quadrupoles) that represent each electronic state and transition densities in the diabatic basis. The DME parameters are obtained through a restrained electrostatic potential (RESP) fit, which we extended to yield higher-order multipoles. We also implemented in SHARC a scheme for achieving roto-translational invariance of LVC models as well as a general quantum mechanics/molecular mechanics (QM/MM) interface, an OpenMM interface, and restraining potentials for simulating liquid droplets. Using thioformaldehyde in water as a test case, we demonstrate that LVC/MM can accurately reproduce the solvation structure and energetics of rigid solutes, with errors on the order of 1-2 kcal/mol compared to a BP86/MM reference. The implementation in SHARC is shown to be very efficient, enabling the simulation of trajectories on the nanosecond time scale in a matter of days.