Force Field Benchmark of Amino Acids: I. Hydration and Diffusion in Different Water Models

Modeling solvation dynamics of transition metal redox ion through on-the-fly multi-objective Bayesian-optimized force field

Modeling solvation dynamics and properties is crucial for developing electrolytes for electrochemical energy storage and conversion devices. This work reports an on-the-fly multi-objective Bayesian optimization (OTF-MOBO) method to parameterize force fields for modeling ionic solvation structures, thermodynamics, and transport properties using molecular dynamics simulations. By leveraging solvation-free energy and solvation radii as training data, we employ the data-driven OTF-MOBO algorithm to actively optimize the force field parameters. The modeling accuracy was evaluated in molecular dynamics simulations until the Pareto front in the parameter space was reached through minimized prediction errors in both solvation-free energy and solvation radii. Using transition metal redox ions (Fe3+/Fe2+, Cr3+/Cr2+, and Cu2+/Cu+) in aqueous solution as examples, we demonstrate that simple force fields combining the Lenard-Jones potential and Coulombic potential can achieve relative error below 2% in both solvation free energy and solvation radii. The optimized force fields can be further extrapolated to predict solvation entropy and diffusivities with relative error below 10% compared with experiments.

Methodological and force field effects in the molecular dynamics-based prediction of binding free energies of host-guest systems

As a contribution to the understanding and rationalization of methodological and modeling effects in recent host-guest SAMPL challenges, using an alchemical molecular dynamics technique we have examined the impact of force field parameterization and ionic strength in connection with guest charge neutralization on computed dissociation free energies in two typical SAMPL heavily charged macrocyclic hosts encapsulating small protonated amines with disparate binding affinities. We have shown that the methodological treatment for host neutralization, with explicit ions or with the background neutralizing plasma in the context of alchemical calculations under periodic boundary conditions, has a moderate effect on the calculated affinities. On the other hand, we have shown that seemingly small differences in the force field parameterization in highly symmetric hosts can produce systematic effects on the structural features that can have a significant impact on the predicted binding affinities.

Modeling Lanthanide Ions in Solution: A Versatile Force Field in Aqueous and Organic Solvents

In this paper, we propose a new nonpolarizable force field for describing the Ln3+ (Ln = lanthanide) series based on a 12-6-4 Lennard-Jones potential. The development of the force field was performed in pure water by adjusting both the ion-oxygen distance and the hydration free energy. This force field accurately reproduces the Ln3+ hydration properties through the series, especially the coordination number that is hardly accessible using a nonpolarizable force field. Then, the validity and the transferability of the current force field were evaluated for two different systems containing Ln3+ in various solvents, namely, 0.1 mol L-1 La(NO3)3 salts in methanol and Eu(NO3)3 salts in solvent organic phases composed of DMDOHEMA molecules in n-heptane. The good agreement between our simulations and the data available in the literature confirms the accuracy of the force field for describing the lanthanide cations in both aqueous and nonaqueous media.

Dealing with Induced Fit, Conformational Selection, and Secondary Poses in Molecular Dynamics Simulations for Reliable Free Energy Predictions

In this study, we have tested the performance of standard molecular dynamics (MD) simulations, replicates of shorter standard MD simulations, and Hamiltonian Replica Exchange (HREM) simulations for the sampling of two macrocyclic hosts for guest delivery, characterized by induced fit (phenyl-based host) and conformation selection (naphthyl-based host) and of the ODR-BRD4(I) drug-receptor system where the ligand can assume two main poses. For the optimization of the HREM simulation, we have proposed and tested an on-the-fly iterative scheme for equalizing the acceptance ratio along the replica progression at a constant replica number resulting in a moderate impact of the sampling efficiency. Concerning standard MD, we have found that, while splitting the total allocated simulation time in short MD replicates can reproduce the sampling efficiency of HREM in the phenyl-based host and in the ODR-BRD4(I) complex, in the naphthyl-based macrocycle, characterized by long-lived metastable states, enhanced sampling techniques are the only viable alternative for a reliable canonical sampling of the rugged conformational landscape.

Assessment of an Electrostatic Energy-Based Charge Model for Modeling the Electrostatic Interactions in Water Solvent

The protein force field based on the restrained electrostatic potential (RESP) charges has limitations in accurately describing hydrogen bonding interactions in proteins. To address this issue, we propose an alternative approach called the electrostatic energy-based charges (EEC) model, which shows improved performance in describing electrostatic interactions (EIs) of hydrogen bonds in proteins. In this study, we further investigate the performance of the EEC model in modeling EIs in water solvent. Our findings demonstrate that the fixed EEC model can effectively reproduce the quantum mechanics/molecular mechanics (QM/MM)-calculated EIs between a water molecule and various water solvent environments. However, to achieve the same level of computational accuracy, the electrostatic potential (ESP) charge model needs to fluctuate according to the electrostatic environment. Our analysis indicates that the requirement for charge adjustments depends on the specific mathematical and physical representation of EIs as a function of the environment for deriving charges. By comparing with widely used empirical water models calibrated to reproduce experimental properties, we confirm that the performance of the EEC model in reproducing QM/MM EIs is similar to that of general purpose TIP4P-like water models such as TIP4P-Ew and TIP4P/2005. When comparing the computed 10,000 distinct EI values within the range of -40 to 0 kcal/mol with the QM/MM results calculated at the MP2/aug-cc-pVQZ/TIP3P level, we noticed that the mean unsigned error (MUE) for the EEC model is merely 0.487 kcal/mol, which is remarkably similar to the MUE values of the TIP4P-Ew (0.63 kcal/mol) and TIP4P/2005 (0.579 kcal/mol) models. However, both the RESP method and the TIP3P model exhibit a tendency to overestimate the EIs, as evidenced by their higher MUE values of 1.761 and 1.293 kcal/mol, respectively. EEC-based molecular dynamics simulations have demonstrated that, when combined with appropriate van der Waals parameters, the EEC model can closely reproduce oxygen-oxygen radial distribution function and density of water, showing a remarkable similarity to the well-established TIP4P-like empirical water models. Our results demonstrate that the EEC model has the potential to build force fields with comparable accuracy to more sophisticated empirical TIP4P-like water models.

Calculation of binding affinity of JAK1 inhibitors via accurately computational estimation

Janus kinase 1 (JAK1) is a tyrosine kinase that is involved in the initiation of responses to a number of different cytokine receptor families. The JAK1-dependent pathway is a therapeutic target, and several JAK inhibitors have been developed thanks to intensive research. However, since the ATP binding sites of JAK family members are quite alike, JAK1 inhibitors can thus be less selective, resulting in unanticipated adverse effects. Despite this, minor variations in the ATP-binding site have been extensively used to find a variety of small compounds with different inhibitory properties. Stronger binding affinity of JAK1 inhibitors is believed to be able to reduce the negative effects, leading to better treatment results. Therefore, a thorough computational search that can effectively identify ligands with extremely high binding affinity for JAK1 to serve as promising inhibitors is required. Here, a method combining steered-molecular dynamic (SMD) simulations with a modified linear interaction energy (LIE) model has been developed to evaluate the binding affinities of known JAK1 inhibitors. The correlation coefficient between the estimated and experimental values was 0.72 and a root-mean-square error was 0.97 kcal•mol-1, revealing that the SMD/LIE method can precisely and quickly predict the binding free energies of JAK1 inhibitors. Furthermore, three marine fungus-derived compounds, namely hansforesters E, hansforesters G and tetroazolemycins B, were identified to be particularly promising JAK1 inhibitors, accordingly. These findings show that the SMD/LIE method has a lot of promise for in silico screening of possible JAK1 inhibitors from a vast number of compounds that are now accessible.Communicated by Ramaswamy H. Sarma.

Screening, characterization and specific binding mechanism of aptamers against human plasminogen Kringle 5

Plasminogen Kringle 5 is one of the most potent cytokines identified to inhibit the proliferation and migration of vascular endothelial cells. Herein, six aptamer candidates that specifically bind to Kringle 5 were generated by the systematic evolution of ligands by exponential enrichment (SELEX). After 10 rounds of screening against Kringle 5, a highly enriched ssDNA pool was sequenced and the representative aptamers were subjected to binding assays to evaluate their affinity and specificity. The preferred aptamer KG-4, which demonstrated a low dissociation constant (K_d) of ∼ 432 nM and excellent selectivity for Kringle 5. A conserved "motif" of eight bases located at the stem-loop intersection, common to the aptamer, was further confirmed as the recognition element for binding with Kringle 5. The bulge formed by the motif and depression on the lysine binding site of Kringle 5 were both located at the binding interface, and the "induced fit" between their structures played a central role in the recognition process. Kringle 5 interacts KG-4 primarily through enthalpy-driven van der Waals forces and hydrogen bond. The key nucleotides A34 and C35 at motif on KG-4 and the positively charged amino acids in the loop 1 and loop 4 regions on Kringle 5 play a major role in the interaction. Furthermore, KG-4 dose-dependently reduced the proliferation inhibition of vascular endothelial cells by Kringle 5 and had a blocking effect on the function of Kringle 5 in inhibiting migration and promoting apoptosis of vascular endothelial cells in vitro. This study put a new light on protein-aptamer binding mechanism and may provide insight into the treatment of ischemic diseases by target depletion of Kringle 5.

Brownian Motion in Optical Tweezers, a Comparison between MD Simulations and Experimental Data in the Ballistic Regime

The four most popular water models in molecular dynamics were studied in large-scale simulations of Brownian motion of colloidal particles in optical tweezers and then compared with experimental measurements in the same time scale. We present the most direct comparison of colloidal polystyrene particle diffusion in molecular dynamics simulations and experimental data on the same time scales in the ballistic regime. The four most popular water models, all of which take into account electrostatic interactions, are tested and compared based on yielded results and resources required. Three different conditions were simulated: a freely moving particle and one in a potential force field with two different strengths based on 1 pN/nm and 10 pN/nm. In all cases, the diameter of the colloidal particle was 50 nm. The acquired data were compared with experimental measurements performed using optical tweezers with position capture rates as high as 125 MHz. The experiments were performed in pure water on polystyrene particles with a 1 μm diameter in special microchannel cells.

Non-covalent interactions from a Quantum Chemical Topology perspective

About half a century after its little-known beginnings, the quantum topological approach called QTAIM has grown into a widespread, but still not mainstream, methodology of interpretational quantum chemistry. Although often confused in textbooks with yet another population analysis, be it perhaps an elegant but somewhat esoteric one, QTAIM has been enriched with about a dozen other research areas sharing its main mathematical language, such as Interacting Quantum Atoms (IQA) or Electron Localisation Function (ELF), to form an overarching approach called Quantum Chemical Topology (QCT). Instead of reviewing the latter's role in understanding non-covalent interactions, we propose a number of ideas emerging from the full consequences of the space-filling nature of topological atoms, and discuss how they (will) impact on interatomic interactions, including non-covalent ones. The architecture of a force field called FFLUX, which is based on these ideas, is outlined. A new method called Relative Energy Gradient (REG) is put forward, which is able, by computation, to detect which fragments of a given molecular assembly govern the energetic behaviour of this whole assembly. This method can offer insight into the typical balance of competing atomic energies both in covalent and non-covalent case studies. A brief discussion on so-called bond critical points is given, highlighting concerns about their meaning, mainly in the arena of non-covalent interactions.

Impacts of targeting different hydration free energy references on the development of ion potentials

Hydration free energy (HFE) as the most important solvation parameter is often targeted in ion model development, even though the reported values differ by dozens of kcal mol^-1 mainly due to the experimentally undetermined HFE of the proton ΔG°(H⁺). The choice of ΔG°(H⁺) obviously affects the hydration of single ions and the relative HFE between the ions with different (magnitude or sign) charges, and the impacts of targeted HFEs on the ion solvation and ion-ion interactions are largely unrevealed. Here we designed point charge models of K⁺, Mg²⁺, Al³⁺, and Cl^- ions targeting a variety of HFE references and then investigated the HFE influences on the simulations of dilute and concentrated ion solutions and of the salt ion pairs in gas, liquid, and solid phases. Targeting one more property of ion-water oxygen distances (IOD) leaves the ion-water binding distance invariant, while the binding strength increases with the decreasing (more negative) HFE of ions as a result of a decrease in ΔG°(H⁺) for the cation and an increase in ΔG°(H⁺) for the anion. The increase in ΔG°(H⁺) leads to strengthened cation-anion interactions and thus to close ion-ion contacts, low osmotic pressures, and small activity derivatives in concentrated ion solutions as well as too stable ion pairs of the salts in different phases. The ion diffusivity and water exchange rates around the ions are simply not HFE dependent but rather more complex. Targeting both the aqueous IOD and salt crystal properties of KCl was also attempted and the comparison between different models indicates the complexity and challenge in obtaining a balanced performance between different phases using classical force fields. Our results also support that a real ΔG°(H⁺) value of -259.8 kcal mol^-1 recommended by Hünenberger and Reif guides ion models to reproduce ion-water and ion-ion interactions reasonably at relatively low salt concentrations. Simulations of a metalloprotein show that a relatively more positive ΔG°(H⁺) for Mg²⁺ model is better for a reasonable description of the metal binding network.

Mechanism-Based Strategy for Optimizing HaloTag Protein Labeling

HaloTag labeling technology has introduced unrivaled potential in protein chemistry and molecular and cellular biology. A wide variety of ligands have been developed to meet the specific needs of diverse applications, but only a single protein tag, DhaAHT, is routinely used for their incorporation. Following a systematic kinetic and computational analysis of different reporters, a tetramethylrhodamine- and three 4-stilbazolium-based fluorescent ligands, we showed that the mechanism of incorporating different ligands depends both on the binding step and the efficiency of the chemical reaction. By studying the different haloalkane dehalogenases DhaA, LinB, and DmmA, we found that the architecture of the access tunnels is critical for the kinetics of both steps and the ligand specificity. We showed that highly efficient labeling with specific ligands is achievable with natural dehalogenases. We propose a simple protocol for selecting the optimal protein tag for a specific ligand from the wide pool of available enzymes with diverse access tunnel architectures. The application of this protocol eliminates the need for expensive and laborious protein engineering.

Molecular simulations of complex carbohydrates and glycoconjugates

Complex carbohydrates (glycans) are the most abundant and versatile biopolymers in nature. The broad diversity of biochemical functions that carbohydrates cover is a direct consequence of the variety of 3D architectures they can adopt, displaying branched or linear arrangements, widely ranging in sizes, and with the highest diversity of building blocks of any other natural biopolymer. Despite this unparalleled complexity, a common denominator can be found in the glycans' inherent flexibility, which hinders experimental characterization, but that can be addressed by high-performance computing (HPC)-based molecular simulations. In this short review, I present and discuss the state-of-the-art of molecular simulations of complex carbohydrates and glycoconjugates, highlighting methodological strengths and weaknesses, important insights through emblematic case studies, and suggesting perspectives for future developments.

501Y.V2 spike protein resists the neutralizing antibody in atomistic simulations

SARS-CoV-2 outbreaks worldwide caused COVID-19 pandemic, which is related to several million deaths. In particular, SARS-CoV-2 Spike (S) protein is a major biological target for COVID-19 vaccine design. Unfortunately, recent reports indicated that Spike (S) protein mutations can lead to antibody resistance. However, understanding the process is limited, especially at the atomic scale. The structural change of S protein and neutralizing antibody fragment (FAb) complexes was thus probed using molecular dynamics (MD) simulations. In particular, the backbone RMSD of the 501Y.V2 complex was significantly larger than that of the wild-type one implying a large structural change of the mutation system. Moreover, the mean of Rg, CCS, and SASA are almost the same when compared two complexes, but the distributions of these values are absolutely different. Furthermore, the free energy landscape of the complexes was significantly changed when the 501Y.V2 variant was induced. The binding pose between S protein and FAb was thus altered. The FAb-binding affinity to S protein was thus reduced due to revealing over steered-MD (SMD) simulations. The observation is in good agreement with the respective experiment that the 501Y.V2 SARS-CoV-2 variant can escape from neutralizing antibody (NAb).

Umbrella Sampling-Based Method to Compute Ligand-Binding Affinity

Many proteins have a solvent-exposed binding cleft, which permits their inhibitors to bind and unbind without significant protein conformation transforms. The binding/unbinding pathways of these protein-inhibitor complexes can be rather straightforwardly sampled by using umbrella sampling (US) simulation methods. During a US simulation, the C_α atoms of the protein are restrained via a harmonic force. The potential of mean force (PMF) along the binding pathway can be estimated by using the weighted histogram analysis method (WHAM). The binding affinity is then computed as the difference in PMF between the binding and unbinding states.

Unbinding ligands from SARS-CoV-2 Mpro via umbrella sampling simulations

The umbrella sampling (US) simulation is demonstrated to be an efficient approach for determining the unbinding pathway and binding affinity to the SARS-CoV-2 Mpro of small molecule inhibitors. The accuracy of US is in the same range as the linear interaction energy (LIE) and fast pulling of ligand (FPL) methods. In detail, the correlation coefficient between US and experiments does not differ from FPL and is slightly smaller than LIE. The root mean square error of US simulations is smaller than that of LIE. Moreover, US is better than FPL and poorer than LIE in classifying SARS-CoV-2 Mpro inhibitors owing to the reciever operating characteristic-area under the curve analysis. Furthermore, the US simulations also provide detailed insights on unbinding pathways of ligands from the binding cleft of SARS-CoV-2 Mpro. The residues Cys44, Thr45, Ser46, Leu141, Asn142, Gly143, Glu166, Leu167, Pro168, Ala191, Gln192 and Ala193 probably play an important role in the ligand dissociation. Therefore, substitutions at these points may change the mechanism of binding of inhibitors to SARS-CoV-2 Mpro.

Binding Energy and Free Energy of Calcium Ion to Calmodulin EF-Hands with the Drude Polarizable Force Field

Calcium ions are important messenger molecules in cells, which bind calcium-binding proteins to trigger many biochemical processes. We constructed four model systems, each containing one EF-hand loop of calmodulin with one calcium ion bound, and investigated the binding energy and free energy of Ca²⁺ by the quantum mechanics symmetry-adapted perturbation theory (SAPT) method and the molecular mechanics with the additive CHARMM36m (C36m) and the polarizable Drude force fields (FFs). Our results show that the explicit introduction of polarizability in the Drude not only yields considerably improved agreement with the binding energy calculated from the SAPT method but is also able to capture each component of the binding energies including electrostatic, induction, exchange, and dispersion terms. However, binding free energies computed with the Drude and the C36m FFs both deviated significantly from the experimental measurements. Detailed analysis indicated that one of main reasons might be that the strong interactions between Ca²⁺ and the side chain nitrogen of Asn/Gln in the Drude FF caused the distorted coordination geometries of calcium. Our work illustrated the importance of polarization in modeling ion-protein interactions and the difficulty in generating accurate and balanced FF models to represent the polarization effects.

Rational Design of Nonbonded Point Charge Models for Monovalent Ions with Lennard-Jones 12-6 Potential

Ions are of central importance in nature, and a variety of potential models was proposed to model ions in different phases for an in-depth exploration of ion-related systems. Here, we developed point charge models of 14 monovalent ions with the traditional 12-6 Lennard-Jones (LJ) potential for use in conjunction with 11 water models of TIP3P, OPC3, SPC/E, SPC/E_b, TIP3P-FB, a99SB-disp, TIP4P-Ew, OPC, TIP4P/2005, TIP4P-D, and TIP4P-FB. The designed models reproduced the real hydration free energy (HFE) of ions and the ion-oxygen distance (IOD) in the first hydration shell accurately and simultaneously, a performance similar to the previously reported 12-6-4 LJ-type ion models (12-6 LJ plus an attractive C₄ term for cations or a repulsive one for anions). This work, along with our previous work on di-, tri-, and tetravalent metal cations (J. Chem. Inf. Model. 2021, 61, 4031-4044; J. Chem. Inf. Model. 2021, 61, 4613-4629), demonstrates the feasibility of the simple 12-6 LJ potential in ion modeling. In order for the 12-6 LJ potential to reproduce both the HFE and IOD, the LJ R parameters need to be close to Shannon's ionic radii for the highly charged cations and to the Stokes's van der Waals (vdW) radii for the monovalent ions. With an additional C₄ term, the R parameters of 12-6-4 LJ ion models agree well with the Stokes's vdW radii and have a more physical meaning. It appears that the C₄ term can be merged into the 12-6 LJ potential by a rational tuning of R and the LJ well depth. Simulations of the osmotic coefficients of alkali chloride solutions and the properties of gaseous and solid alkali halides indicate the necessity of further optimizing ion-ion interactions via, for instance, targeting more properties or using a more physical (polarizable) model.

Rational Design of Nonbonded Point Charge Models for Highly Charged Metal Cations with Lennard-Jones 12-6 Potential

Here, we developed nonbonded point charge models using a simple Lennard-Jones (LJ) 12-6 potential for highly charged metal cations (18 trivalent and 6 tetravalent ions) for use with 11 water models of TIP3P, OPC3, SPC/E, SPC/E_b, TIP3P-FB, a99SB-disp, TIP4P-Ew, OPC, TIP4P/2005, TIP4P-D, and TIP4P-FB. The designed models simultaneously reproduce the hydration free energy (HFE) and ion-oxygen distance (IOD) in the first hydration shell with an error within 1 kcal/mol and 0.01 Å on average, respectively, and yield reasonable coordination numbers for most cations. Such performance is equivalent to the previously reported point charge models using a more complex 12-6-4 LJ-type potential, while the LJ R parameters of our models are much close to Shannon's revised effective ion radii than that of the 12-6-4 models. Our designed models overestimate the diffusion constants of several trivalent ions by 5-68%. The performance in predicting osmotic coefficients of trivalent chlorides in aqueous solution depends on the salt type. A calibration of cation-anion interacting LJ parameters reproduces the experimental osmotic coefficients of an AlCl₃ solution at 0.2-3.0 mol/L. The effectiveness of our new models is further demonstrated by simulating a metalloprotein system with four force field/water combinations. This work facilitates accurate modeling of metal-containing systems by a variety of force fields and water models in aqueous solution.

Systematic Comparison of the Structural and Dynamic Properties of Commonly Used Water Models for Molecular Dynamics Simulations

Water is a unique solvent that is ubiquitous in biology and present in a variety of solutions, mixtures, and materials settings. It therefore forms the basis for all molecular dynamics simulations of biological phenomena, as well as for many chemical, industrial, and materials investigations. Over the years, many water models have been developed, and it remains a challenge to find a single water model that accurately reproduces all experimental properties of water simultaneously. Here, we report a comprehensive comparison of structural and dynamic properties of 30 commonly used 3-point, 4-point, 5-point, and polarizable water models simulated using consistent settings and analysis methods. For the properties of density, coordination number, surface tension, dielectric constant, self-diffusion coefficient, and solvation free energy of methane, models published within the past two decades consistently show better agreement with experimental values compared to models published earlier, albeit with some notable exceptions. However, no single model reproduced all experimental values exactly, highlighting the need to carefully choose a water model for a particular study, depending on the phenomena of interest. Finally, machine learning algorithms quantified the relationship between the water model force field parameters and the resulting bulk properties, providing insight into the parameter-property relationship and illustrating the challenges of developing a water model that can accurately reproduce all properties of water simultaneously.

Emergence of Electric Fields at the Water-C12E6 Surfactant Interface

We study the properties of the interface of water and the surfactant hexaethylene glycol monododecyl ether (C12E6) with a combination of heterodyne-detected vibrational sum frequency generation (HD-VSFG), Kelvin-probe measurements, and molecular dynamics (MD) simulations. We observe that the addition of the hydrogen-bonding surfactant C12E6, close to the critical micelle concentration (CMC), induces a drastic enhancement in the hydrogen bond strength of the water molecules close to the interface, as well as a flip in their net orientation. The mutual orientation of the water and C12E6 molecules leads to the emergence of a broad (∼3 nm) interface with a large electric field of ∼1 V/nm, as evidenced by the Kelvin-probe measurements and MD simulations. Our findings may open the door for the design of novel electric-field-tuned catalytic and light-harvesting systems anchored at the water-surfactant-air interface.

Marine derivatives prevent wMUS81 in silico studies

The winged-helix domain of the methyl methanesulfonate and ultraviolet-sensitive 81 (wMUS81) is a potential cancer drug target. In this context, marine fungi compounds were indicated to be able to prevent wMUS81 structure via atomistic simulations. Eight compounds such as D197 (Tryptoquivaline U), D220 (Epiremisporine B), D67 (Aspergiolide A), D153 (Preussomerin G), D547 (12,13-dihydroxyfumitremorgin C), D152 (Preussomerin K), D20 (Marinopyrrole B) and D559 (Fumuquinazoline K) were indicated that they are able to prevent the conformation of wMUS81 via forming a strong binding affinity to the enzyme via perturbation approach. The electrostatic interaction is the dominant factor in the binding process of ligands to wMUS81. The residues Trp55, Arg59, Leu62, His63 and Arg69 were found to frequently form non-bonded contacts and hydrogen bonds to inhibitors. Moreover, the influence of the ligand D197, which formed the lowest binding free energy to wMUS81, on the structural change of enzyme was investigated using replica exchange molecular dynamics simulations. The obtained results indicated that D197, which forms a strong binding affinity, can modify the structure of wMUS81. Overall, the marine compounds probably inhibit wMUS81 due to forming a strong binding affinity to the enzyme as well as altering the enzymic conformation.

Rational Design of Nonbonded Point Charge Models for Divalent Metal Cations with Lennard-Jones 12-6 Potential

Exploring a metal-involved biochemical process at a molecular level often requires a reliable description of metal properties in aqueous solution by classical nonbonded models. An additional C₄ term for considering ion-induced dipole interactions was previously proposed to supplement the widely used Lennard-Jones 12-6 potential (known as the 12-6-4 LJ-type model) with good accuracy. Here, we demonstrate an alternative to modeling divalent metal cations (M²⁺) with the traditional 12-6 LJ potential by developing nonbonded point charge models for use with 11 water models: TIP3P, SPC/E, SPC/E_b, TIP4P-Ew, TIP4P-D, and TIP4P/2005 and the more recent OPC3, TIP3P-FB, OPC, TIP4P-FB, and a99SB-disp. Our designed models simultaneously reproduce the experimental hydration free energy, ion-oxygen distance, and coordination number in the first hydration shell accurately for most of the metal cations, an accuracy equivalent to that of the complex 12-6-4 LJ-type and double exponential potential models. A systematic comparison with the existing M²⁺ models is presented as well in terms of effective ion radii, diffusion constants, water exchange rates, and ion-water interactions. Molecular dynamics simulations of metal substitution in Escherichia coli glyoxalase I variants show the great potential of our new models for metalloproteins.

Force Field Benchmark of Amino Acids. 3. Hydration with Scaled Lennard-Jones Interactions

Classical protein force fields were reported with too weak protein-water interactions relative to protein-protein interactions, leading to more compact structures and artificial protein aggregation. Here we investigated the impacts of scaled Lennard-Jones (LJ) interactions on the hydration of amino acids and the simulation of folded and intrinsically disordered proteins (IDPs). The obtained optimal scaling parameters reproduce accurately hydration free energies of neutral amino acid side chain analogues and do not affect the compactness and structural stability of folded proteins significantly. The scaling leads to less compact IDPs and varies from case to case. Strengthening the interactions between protein and water oxygen or hydrogen atoms by increasing the interacting LJ well depth (ε) appears more effective than weakening protein-protein interactions by reducing the interacting dispersion coefficients (C₆). We demonstrate that weakening water-water interactions is a solution as well to obtaining more favorable protein-water interactions in an indirect way, although modern force fields like Amber ff19SB and a99SB-disp tend to use water models with strong water-water interactions. This is likely a compromise between strong protein-protein interactions and strong water-water interactions. Independent optimization of protein force fields and water models is therefore needed to make both interactions more close to reality, leading to good accuracy without bias or scaling.

Accurate Structure Prediction for Protein Loops Based on Molecular Dynamics Simulations with RSFF2C

Protein loops, connecting the α-helices and β-strands, are involved in many important biological processes. However, due to their conformational flexibility, it is still challenging to accurately determine three-dimensional (3D) structures of long loops experimentally and computationally. Herein, we present a systematic study of the protein loop structure prediction via a total of ∼850 μs molecular dynamics (MD) simulations. For a set of 15 long (10-16 residues) and solvent-exposed loops, we first evaluated the performance of four state-of-the-art loop modeling algorithms, DaReUS-Loop, Sphinx, Rosetta-NGK, and MODELLER, on each loop, and none of them could accurately predict the structures for most loops. Then, temperature replica exchange molecular dynamics (REMD) simulations were conducted with three recent force fields, RSFF2C with TIP3P water model, CHARMM36m with CHARMM-modified TIP3P, and AMBER ff19SB with OPC. We found that our recently developed residue-specific force field RSFF2C performed the best and successfully predicted 12 out of 15 loops with a root-mean-square deviation (RMSD) < 1.5 Å. As an alternative with lower computational cost, normal MD simulations at high temperatures (380, 500, and 620 K) were investigated. Temperature-dependent performance was observed for each force field, and, for RSFF2C+TIP3P, we found that three independent 100-ns MD simulations at 500 K gave comparable results with REMD simulations. These results suggest that MD simulations, especially with enhanced sampling techniques such as replica exchange, with the RSFF2C force field could be useful for accurate loop structure prediction.

Evaluation of nine condensed-phase force fields of the GROMOS, CHARMM, OPLS, AMBER, and OpenFF families against experimental cross-solvation free energies

Experimental solvation free energies are nowadays commonly included as target properties in the validation of condensed-phase force fields, sometimes even in their calibration. In a previous article [Kashefolgheta et al., J. Chem. Theory. Comput., 2020, 16, 7556-7580], we showed how the involved comparison between experimental and simulation results could be made more systematic by considering a full matrix of cross-solvation free energies . For a set of N molecules that are all in the liquid state under ambient conditions, such a matrix encompasses N×N entries for considering each of the N molecules either as solute (A) or as solvent (B). In the quoted study, a cross-solvation matrix of 25 × 25 experimental value was introduced, considering 25 small molecules representative for alkanes, chloroalkanes, ethers, ketones, esters, alcohols, amines, and amides. This experimental data was used to compare the relative accuracies of four popular condensed-phase force fields, namely GROMOS-2016H66, OPLS-AA, AMBER-GAFF, and CHARMM-CGenFF. In the present work, the comparison is extended to five additional force fields, namely GROMOS-54A7, GROMOS-ATB, OPLS-LBCC, AMBER-GAFF2, and OpenFF. Considering these nine force fields, the correlation coefficients between experimental values and simulation results range from 0.76 to 0.88, the root-mean-square errors (RMSEs) from 2.9 to 4.8 kJ mol-1, and average errors (AVEEs) from -1.5 to +1.0 kJ mol-1. In terms of RMSEs, GROMOS-2016H66 and OPLS-AA present the best accuracy (2.9 kJ mol-1), followed by OPLS-LBCC, AMBER-GAFF2, AMBER-GAFF, and OpenFF (3.3 to 3.6 kJ mol-1), and then by GROMOS-54A7, CHARM-CGenFF, and GROMOS-ATB (4.0 to 4.8 kJ mol-1). These differences are statistically significant but not very pronounced, and are distributed rather heterogeneously over the set of compounds within the different force fields.

Molecular dynamics simulation for membrane separation and porous materials: A current state of art review

Computational frameworks have been under specific attention within the last two decades. Molecular Dynamics (MD) simulations, identical to the other computational approaches, try to address the unknown question, lighten the dark areas of unanswered questions, to achieve probable explanations and solutions. Owing to their complex microporous structure on one side and the intricate biochemical nature of various materials used in the structure, separative membrane materials possess peculiar degrees of complications. More notably, as nanocomposite materials are often integrated into separative membranes, thin-film nanocomposites and porous separative nanocomposite materials could possess an additional level of complexity with regard to the nanoscale interactions brought to the structure. This critical review intends to cover the recent methods used to assess membranes and membrane materials. Incorporation of MD in membrane technology-related fields such as desalination, fuel cell-based energy production, blood purification through hemodialysis, etc., were briefly covered. Accordingly, this review could be used to understand the current extent of MD applications for separative membranes. The review could also be used as a guideline to use the proper MD implementation within the related fields.

Molecular design of anticancer drugs from marine fungi derivatives

Heat shock protein 90 (Hsp90) is one of the most potential targets in cancer therapy. We have demonstrated using a combination of molecular docking and fast pulling of ligand (FPL) simulations that marine fungi derivatives can be possible inhibitors, preventing the biological activity of Hsp90. The computational approaches were validated and compared with previous experiments. Based on the benchmark of available inhibitors of Hsp90, the GOLD docking package using the ChemPLP scoring function was found to be superior over both Autodock Vina and Autodock4 in the preliminary estimation of the ligand-binding affinity and binding pose with the Pearson correlation, R = -0.62. Moreover, FPL calculations were also indicated as a suitable approach to refine docking simulations with a correlation coefficient with the experimental data of R = -0.81. Therefore, the binding affinity of marine fungi derivatives to Hsp90 was evaluated. Docking and FPL calculations suggest that five compounds including 23, 40, 46, 48, and 52 are highly potent inhibitors for Hsp90. The obtained results enhance cancer therapy research.

Developing and Assessing Nonbonded Dummy Models of Magnesium Ion with Different Hydration Free Energy References

A large diversity in the targeted hydration free energies (HFEs) during model parameterization of metal ions was reported in the literature with a difference by dozens of kcal/mol. Here, we developed a series of nonbonded dummy models of the Mg²⁺ ion targeting different HFE references in TIP3P water, followed by assessments of the designed models in the simulations of MgCl₂ solution and biological systems. Together with the comparison of existing models, we conclude that the difference in the targeted HFEs has a limited influence on the model performance, while the usability of these models differs from case to case. The feasibility of reproducing more properties of Mg²⁺ such as diffusion constants and water exchange rates using a nonbonded dummy model is demonstrated. Underestimated activity derivative and osmotic coefficient of MgCl₂ solutions in high concentration reveal a necessity for further optimization of ion-pair interactions. The developed dummy models are applicable to metal coordination with Asp, Glu, and His residues in metalloenzymes, and the performance in predicting monodentate or bidentate binding modes of Asp/Glu residues depends on the complexity of metal centers and the choice of protein force fields. When both the binding modes coexist, the nonbonded dummy models outperform point charge models, probably in need of considering polarization of metal-binding residues by, for instance, charge calibration in classical force fields. This work is valuable for the use and further development of magnesium ion models for simulations of metal-containing systems with good accuracy.

Proteins maintain hydration at high [KCl] concentration regardless of content in acidic amino acids

Proteins of halophilic organisms, which accumulate molar concentrations of KCl in their cytoplasm, have a much higher content in acidic amino acids than proteins of mesophilic organisms. It has been proposed that this excess is necessary to maintain proteins hydrated in an environment with low water activity, either via direct interactions between water and the carboxylate groups of acidic amino acids or via cooperative interactions between acidic amino acids and hydrated cations. Our simulation study of five halophilic proteins and five mesophilic counterparts does not support either possibility. The simulations use the AMBER ff14SB force field with newly optimized Lennard-Jones parameters for the interactions between carboxylate groups and potassium ions. We find that proteins with a larger fraction of acidic amino acids indeed have higher hydration levels, as measured by the concentration of water in their hydration shell and the number of water/protein hydrogen bonds. However, the hydration level of each protein is identical at low (b_KCl = 0.15 mol/kg) and high (b_KCl = 2 mol/kg) KCl concentrations; excess acidic amino acids are clearly not necessary to maintain proteins hydrated at high salt concentration. It has also been proposed that cooperative interactions between acidic amino acids in halophilic proteins and hydrated cations stabilize the folded protein structure and would lead to slower dynamics of the solvation shell. We find that the translational dynamics of the solvation shell is barely distinguishable between halophilic and mesophilic proteins; if such a cooperative effect exists, it does not have that entropic signature.

Benchmark of Popular Free Energy Approaches Revealing the Inhibitors Binding to SARS-CoV-2 Mpro

The COVID-19 pandemic has killed millions of people worldwide since its outbreak in December 2019. The pandemic is caused by the SARS-CoV-2 virus whose main protease (Mpro) is a promising drug target since it plays a key role in viral proliferation and replication. Currently, developing an effective therapy is an urgent task, which requires accurately estimating the ligand-binding free energy to SARS-CoV-2 Mpro. However, it should be noted that the accuracy of a free energy method probably depends on the protein target. A highly accurate approach for some targets may fail to produce a reasonable correlation with the experiment when a novel enzyme is considered as a drug target. Therefore, in this context, the ligand-binding affinity to SARS-CoV-2 Mpro was calculated via various approaches. The molecular docking approach was manipulated using Autodock Vina (Vina) and Autodock4 (AD4) protocols to preliminarily investigate the ligand-binding affinity and pose to SARS-CoV-2 Mpro. The binding free energy was then refined using the fast pulling of ligand (FPL), linear interaction energy (LIE), molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA), and free energy perturbation (FEP) methods. The benchmark results indicated that for docking calculations, Vina is more accurate than AD4, and for free energy methods, FEP is the most accurate method, followed by LIE, FPL, and MM-PBSA (FEP > LIE ≈ FPL > MM-PBSA). Moreover, atomistic simulations revealed that the van der Waals interaction is the dominant factor. The residues Thr26, His41, Ser46, Asn142, Gly143, Cys145, His164, Glu166, and Gln189 are essential elements affecting the binding process. Our benchmark provides guidelines for further investigations using computational approaches.

Computational estimation of potential inhibitors from known drugs against the main protease of SARS-CoV-2

The coronavirus disease (COVID-19) pandemic caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has rapidly spread worldwide recently, leading to global social and economic disruption. Although the emergently approved vaccine programs against SARS-CoV-2 have been rolled out globally, the number of COVID-19 daily cases and deaths has remained significantly high. Here, we attempt to computationally screen for possible medications for COVID-19 via rapidly estimating the highly potential inhibitors from an FDA-approved drug database against the main protease (Mpro) of SARS-CoV-2. The approach combined molecular docking and fast pulling of ligand (FPL) simulations that were demonstrated to be accurate and suitable for quick prediction of SARS-CoV-2 Mpro inhibitors. The results suggested that twenty-seven compounds were capable of strongly associating with SARS-CoV-2 Mpro. Among them, the seven top leads are daclatasvir, teniposide, etoposide, levoleucovorin, naldemedine, cabozantinib, and irinotecan. The potential application of these drugs in COVID-19 therapy has thus been discussed.

Not only in silico drug discovery: Molecular modeling towards in silico drug delivery formulations

The use of methods at molecular scale for the discovery of new potential active ligands, as well as previously unknown binding sites for target proteins, is now an established reality. Literature offers many successful stories of active compounds developed starting from insights obtained in silico and approved by Food and Drug Administration (FDA). One of the most famous examples is raltegravir, a HIV integrase inhibitor, which was developed after the discovery of a previously unknown transient binding area thanks to molecular dynamics simulations. Molecular simulations have the potential to also improve the design and engineering of drug delivery devices, which are still largely based on fundamental conservation equations. Although they can highlight the dominant release mechanism and quantitatively link the release rate to design parameters (size, drug loading, et cetera), their spatial resolution does not allow to fully capture how phenomena at molecular scale influence system behavior. In this scenario, the "computational microscope" offered by simulations at atomic scale can shed light on the impact of molecular interactions on crucial parameters such as release rate and the response of the drug delivery device to external stimuli, providing insights that are difficult or impossible to obtain experimentally. Moreover, the new paradigm brought by nanomedicine further underlined the importance of such computational microscope to study the interactions between nanoparticles and biological components with an unprecedented level of detail. Such knowledge is a fundamental pillar to perform device engineering and to achieve efficient and safe formulations. After a brief theoretical background, this review aims at discussing the potential of molecular simulations for the rational design of drug delivery systems.

Characterizing the Water Wire in the Gramicidin Channel Found by Monte Carlo Sampling Using Continuum Electrostatics and in Molecular Dynamics Trajectories with Conventional or Polarizable Force Fields

Water molecules play a key role in all biochemical processes. They help define the shape of proteins, and they are reactant or product in many reactions and are released as ligands are bound. They facilitate the transfer of protons through transmembrane proton channel, pump and transporter proteins. Continuum electrostatics (CE) force fields used by program Multiconformation CE (MCCE) capture electrostatic interactions in biomolecules with an implicit solvent, which captures the averaged solvent water equilibrium properties. Hybrid CE methods can use explicit water molecules within the protein surrounded by implicit solvent. These hybrid methods permit the study of explicit hydrogen bond networks within the protein and allow analysis of processes such as proton transfer reactions. Yet hybrid CE methods have not been rigorously tested. Here, we present an explicit treatment of water molecules in the Gramicidin A (gA) channel using MCCE and compare the resulting distributions of water molecules and key hydration features against those obtained with explicit solvent Molecular Dynamics (MD) simulations with the nonpolarizable CHARMM36 and polarizable Drude force fields. CHARMM36 leads to an aligned water wire in the channel characterized by a large absolute net water dipole moment; the MCCE and Drude analysis lead to a small net dipole moment as the water molecules change orientation within the channel. The correct orientation is not as yet known, so these calculations identify an open question.

Molecular Simulations of Hydrophobic Gating of Pentameric Ligand Gated Ion Channels: Insights into Water and Ions

Ion channels are proteins which form gated nanopores in biological membranes. Many channels exhibit hydrophobic gating, whereby functional closure of a pore occurs by local dewetting. The pentameric ligand gated ion channels (pLGICs) provide a biologically important example of hydrophobic gating. Molecular simulation studies comparing additive vs polarizable models indicate predictions of hydrophobic gating are robust to the model employed. However, polarizable models suggest favorable interactions of hydrophobic pore-lining regions with chloride ions, of relevance to both synthetic carriers and channel proteins. Electrowetting of a closed pLGIC hydrophobic gate requires too high a voltage to occur physiologically but may inform designs for switchable nanopores. Global analysis of ∼200 channels yields a simple heuristic for structure-based prediction of (closed) hydrophobic gates. Simulation-based analysis is shown to provide an aid to interpretation of functional states of new channel structures. These studies indicate the importance of understanding the behavior of water and ions within the nanoconfined environment presented by ion channels.

Computational design of enzymes for biotechnological applications

Enzymes are the natural catalysts that execute biochemical reactions upholding life. Their natural effectiveness has been fine-tuned as a result of millions of years of natural evolution. Such catalytic effectiveness has prompted the use of biocatalysts from multiple sources on different applications, including the industrial production of goods (food and beverages, detergents, textile, and pharmaceutics), environmental protection, and biomedical applications. Natural enzymes often need to be improved by protein engineering to optimize their function in non-native environments. Recent technological advances have greatly facilitated this process by providing the experimental approaches of directed evolution or by enabling computer-assisted applications. Directed evolution mimics the natural selection process in a highly accelerated fashion at the expense of arduous laboratory work and economic resources. Theoretical methods provide predictions and represent an attractive complement to such experiments by waiving their inherent costs. Computational techniques can be used to engineer enzymatic reactivity, substrate specificity and ligand binding, access pathways and ligand transport, and global properties like protein stability, solubility, and flexibility. Theoretical approaches can also identify hotspots on the protein sequence for mutagenesis and predict suitable alternatives for selected positions with expected outcomes. This review covers the latest advances in computational methods for enzyme engineering and presents many successful case studies.

Ligand Affinities within the Open-Boundary Molecular Mechanics/Coarse-Grained Framework (I): Alchemical Transformations within the Hamiltonian Adaptive Resolution Scheme

Our recently developed Open-Boundary Molecular Mechanics/Coarse Grained (OB-MM/CG) framework predicts ligand poses in important pharmaceutical targets, such as G-protein Coupled Receptors, even when experimental structural information is lacking. The approach, which is based on GROMOS and AMBER force fields, allows for grand-canonical simulations of protein-ligand complexes by using the Hamiltonian Adaptive Resolution Scheme (H-AdResS) for the solvent. Here, we present a key step toward the estimation of ligand binding affinities for their targets within this approach. This is the implementation of the H-AdResS in the GROMACS code. The accuracy of our implementation is established by calculating hydration free energies of several molecules in water by means of alchemical transformations. The deviations of the GROMOS- and AMBER-based H-AdResS results from the reference fully atomistic simulations are smaller than the accuracy of the force field and/or they are in the range of the published results. Importantly, our predictions are in good agreement with experimental data. The current implementation paves the way to the use of the OB-MM/CG framework for the study of large biological systems.

Methodological uncertainties in drug-receptor binding free energy predictions based on classical molecular dynamics

Computational approaches are becoming an essential tool in modern drug design and discovery, with fast compound triaging using a combination of machine learning and docking techniques followed by molecular dynamics binding free energies assessment using alchemical techniques. The traditional MD-based alchemical free energy perturbation (FEP) method faces severe sampling issues that may limits its reliability in automated workflows. Here we review the major sources of uncertainty in FEP protocols for drug discovery, showing how the sampling problem can be effectively tackled by switching to nonequilibrium alchemical techniques.

Estimating the ligand-binding affinity via λ-dependent umbrella sampling simulations

The umbrella sampling (US) approach has been demonstrated to be a very efficient method for estimating the ligand-binding affinity. However, most of the calculated values overestimate experimental ones that are probably caused by the inaccurate representation of the interaction between the ligand and the surrounding molecules. The issue can be resolved via the implementation aspects of λ-alteration simulation into the US approach, which we call the λ-dependent umbrella sampling (λUS) scheme. In particular, the electrostatic and van der Waals interactions were simultaneously changed by using the coupling parameter λ during λUS simulations. The mean value of obtained results, ∆ G US λ = 0.20 = - 11.59 ± 1.51 kcal mol^-1 , is in good fitting to the mean value of respective experiments, ∆G_EXP =  - 11.26 ± 0.89 kcal mol^-1 . Moreover, the correlation between the proposed approach and experiment is quite good with a value of R US λ = 0.20 = 0.82 ± 0.10 . The λUS scheme significantly enhances the calculated accuracy since the RMSE of the proposed scheme is smaller than traditional US simulations, RMSE US λ = 0.20 = 2.99 ± 0.82 kcal mol^-1 versus RMSE US λ = 0.00 = 5.48 ± 0.81 kcal mol^-1 . Furthermore, the precision is increased since the computed error via λUS approach, δ US λ = 0.20 = 1.51 kcal mol^-1 , was smaller than those of the US simulation, δ US λ = 0.00 = 1.78 kcal mol^-1 . Overall, the proposed approach perhaps provides an efficient way to accurately and precisely estimate the ligand-binding free energy.

Systematic design of biomolecular force fields

Force fields for the study of biomolecules have been developed in a predominantly organic manner by regular updates over half a century. Together with better algorithms and advances in computer technology, force fields have improved to yield more robust predictions. However, there are also indications to suggest that intramolecular energy functions have not become better and that there still is room for improvement. In this review, systematic efforts toward development of novel force fields from scratch are described. This includes an estimate of the complexity of the problem and the prerequisites in the form of data and algorithms. It is suggested that in order to make progress, an effort is needed to standardize reference data and force field validation benchmarks.

Water in Nanopores and Biological Channels: A Molecular Simulation Perspective

This Review explores the dynamic behavior of water within nanopores and biological channels in lipid bilayer membranes. We focus on molecular simulation studies, alongside selected structural and other experimental investigations. Structures of biological nanopores and channels are reviewed, emphasizing those high-resolution crystal structures, which reveal water molecules within the transmembrane pores, which can be used to aid the interpretation of simulation studies. Different levels of molecular simulations of water within nanopores are described, with a focus on molecular dynamics (MD). In particular, models of water for MD simulations are discussed in detail to provide an evaluation of their use in simulations of water in nanopores. Simulation studies of the behavior of water in idealized models of nanopores have revealed aspects of the organization and dynamics of nanoconfined water, including wetting/dewetting in narrow hydrophobic nanopores. A survey of simulation studies in a range of nonbiological nanopores is presented, including carbon nanotubes, synthetic nanopores, model peptide nanopores, track-etched nanopores in polymer membranes, and hydroxylated and functionalized nanoporous silica. These reveal a complex relationship between pore size/geometry, the nature of the pore lining, and rates of water transport. Wider nanopores with hydrophobic linings favor water flow whereas narrower hydrophobic pores may show dewetting. Simulation studies over the past decade of the behavior of water in a range of biological nanopores are described, including porins and β-barrel protein nanopores, aquaporins and related polar solute pores, and a number of different classes of ion channels. Water is shown to play a key role in proton transport in biological channels and in hydrophobic gating of ion channels. An overall picture emerges, whereby the behavior of water in a nanopore may be predicted as a function of its hydrophobicity and radius. This informs our understanding of the functions of diverse channel structures and will aid the design of novel nanopores. Thus, our current level of understanding allows for the design of a nanopore which promotes wetting over dewetting or vice versa. However, to design a novel nanopore, which enables fast, selective, and gated flow of water de novo would remain challenging, suggesting a need for further detailed simulations alongside experimental evaluation of more complex nanopore systems.

Rapid prediction of possible inhibitors for SARS-CoV-2 main protease using docking and FPL simulations

Originating for the first time in Wuhan, China, the outbreak of SARS-CoV-2 has caused a serious global health issue. An effective treatment for SARS-CoV-2 is still unavailable. Therefore, in this study, we have tried to predict a list of potential inhibitors for SARS-CoV-2 main protease (Mpro) using a combination of molecular docking and fast pulling of ligand (FPL) simulations. The approaches were initially validated over a set of eleven available inhibitors. Both Autodock Vina and FPL calculations produced consistent results with the experiments with correlation coefficients of R Dock = 0.72 ± 0.14 and R W = -0.76 ± 0.10, respectively. The combined approaches were then utilized to predict possible inhibitors that were selected from a ZINC15 sub-database for SARS-CoV-2 Mpro. Twenty compounds were suggested to be able to bind well to SARS-CoV-2 Mpro. Among them, five top-leads are periandrin V, penimocycline, cis-p-Coumaroylcorosolic acid, glycyrrhizin, and uralsaponin B. The obtained results could probably lead to enhance the COVID-19 therapy.

Estimation of the ligand-binding free energy of checkpoint kinase 1 via non-equilibrium MD simulations

Checkpoint kinase 1 (CHK1) is a serine/threonine-protein kinase that is involved in cell cycle regulation in eukaryotes. Inhibition of CHK1 is thus considered as a promising approach in cancer therapy. In this study, the fast pulling of ligand (FPL) process was applied to predict the relative binding affinities of CHK1 inhibitors using non-equilibrium molecular dynamics (MD) simulations. The work of external harmonic forces to pull the ligand out of the binding cavity strongly correlated with the experimental binding affinity of CHK1 inhibitors with the correlation coefficient of R = -0.88 and an overall root mean square error (RMSE) of 0.99 kcal/mol. The data indicate that the FPL method is highly accurate in predicting the relative binding free energies of CHK1 inhibitors with an affordable CPU time. A new set of molecules were designed based on the molecular modeling of interactions between the known inhibitor and CHK1 as inhibitory candidates. Molecular docking and FPL results exhibited that the binding affinities of developed ligands were similar to the known inhibitor in interaction with the catalytic site of CHK1, producing very potential CHK1 inhibitors of that the inhibitory activities should be further evaluated in vitro.

Induced Polarization in Molecular Dynamics Simulations of the 5-HT3 Receptor Channel

Ion channel proteins form water-filled nanoscale pores within lipid bilayers, and their properties are dependent on the complex behavior of water in a nanoconfined environment. Using a simplified model of the pore of the 5-HT3 receptor (5HT3R) which restrains the backbone structure to that of the parent channel protein from which it is derived, we compare additive with polarizable models in describing the behavior of water in nanopores. Molecular dynamics simulations were performed with four conformations of the channel: two closed state structures, an intermediate state, and an open state, each embedded in a phosphatidylcholine bilayer. Water density profiles revealed that for all water models, the closed and intermediate states exhibited strong dewetting within the central hydrophobic gate region of the pore. However, the open state conformation exhibited varying degrees of hydration, ranging from partial wetting for the TIP4P/2005 water model to complete wetting for the polarizable AMOEBA14 model. Water dipole moments calculated using polarizable force fields also revealed that water molecules remaining within dewetted sections of the pore resemble gas phase water. Free energy profiles for Na+ and for Cl- ions within the open state pore revealed more rugged energy landscapes using polarizable force fields, and the hydration number profiles of these ions were also sensitive to induced polarization resulting in a substantive reduction of the number of waters within the first hydration shell of Cl- while it permeates the pore. These results demonstrate that induced polarization can influence the complex behavior of water and ions within nanoscale pores and provides important new insights into their chemical properties.

Atomic Partitioning of the MPn (n = 2, 3, 4) Dynamic Electron Correlation Energy by the Interacting Quantum Atoms Method: A Fast and Accurate Electrostatic Potential Integral Approach

Recently, the quantum topological energy partitioning method called interacting quantum atoms (IQA) has been extended to MPn (n = 2, 3, 4) wave functions. This enables the extraction of chemical insight related to dynamic electron correlation. The large computational expense of the IQA-MPn approach is compensated by the advantages that IQA offers compared to older nontopological energy decomposition schemes. This expense is problematic in the construction of a machine learning training set to create kriging models for topological atoms. However, the algorithm presented here markedly accelerates the calculation of atomically partitioned electron correlation energies. Then again, the algorithm cannot calculate pairwise interatomic energies because it applies analytical integrals over whole space (rather than over atomic volumes). However, these pairwise energies are not needed in the quantum topological force field FFLUX, which only uses the energy of an atom interacting with all remaining atoms of the system that it is part of. Thus, it is now feasible to generate accurate and sizeable training sets at MPn level of theory. © 2019 The Authors. Journal of Computational Chemistry published by Wiley Periodicals, Inc.

Rotational and Translational Diffusion of Proteins as a Function of Concentration

Atomistic simulations of three different proteins at different concentrations are performed to obtain insight into protein mobility as a function of protein concentration. We report on simulations of proteins from diluted to the physiological water concentration (about 70% of the mass). First, the viscosity was computed and found to increase by a factor of 7-9 going from pure water to the highest protein concentration, in excellent agreement with in vivo nuclear magnetic resonance results. At a physiological concentration of proteins, the translational diffusion is found to be slowed down to about 30% of the in vitro values. The slow-down of diffusion found here using atomistic models is slightly more than that of a hard sphere model that neglects the electrostatic interactions. Interestingly, rotational diffusion of proteins is slowed down somewhat more (by about 80-95% compared to in vitro values) than translational diffusion, in line with experimental findings and consistent with the increased viscosity. The finding that rotation is retarded more than translation is attributed to solvent-separated clustering. No direct interactions between the proteins are found, and the clustering can likely be attributed to dispersion interactions that are stronger between proteins than between protein and water. Based on these simulations, we can also conclude that the internal dynamics of the proteins in our study are affected only marginally under crowding conditions, and the proteins become somewhat more stable at higher concentrations. Simulations were performed using a force field that was tuned for dealing with crowding conditions by strengthening the protein-water interactions. This force field seems to lead to a reproducible partial unfolding of an α-helix in one of the proteins, an effect that was not observed in the unmodified force field.