Parametrization of Trivalent and Tetravalent Metal Ions for the OPC3, OPC, TIP3P-FB, and TIP4P-FB Water Models

Simulating Metal-Imidazole Complexes

One commonly observed binding motif in metalloproteins involves the interaction between a metal ion and histidine's imidazole side chains. Although previous imidazole-M(II) parameters established the flexibility and reliability of the 12-6-4 Lennard-Jones (LJ)-type nonbonded model by simply tuning the ligating atom's polarizability, they have not been applied to multiple-imidazole complexes. To fill this gap, we systematically simulate multiple-imidazole complexes (ranging from one to six) for five metal ions (Co(II), Cu(II), Mn(II), Ni(II), and Zn(II)) which commonly appear in metalloproteins. Using extensive (40 ns per PMF window) sampling to assemble free energy association profiles (using OPC water and standard HID imidazole charge models from AMBER) and comparing the equilibrium distances to DFT calculations, a new set of parameters was developed to focus on energetic and geometric features of multiple-imidazole complexes. The obtained free energy profiles agree with the experimental binding free energy and DFT calculated distances. To validate our model, we show that we can close the thermodynamic cycle for metal-imidazole complexes with up to six imidazole molecules in the first solvation shell. Given the success in closing the thermodynamic cycles, we then used the same extended sampling method for six other metal ions (Ag(I), Ca(II), Cd(II), Cu(I), Fe(II), and Mg(II)) to obtain new parameters. Since these new parameters can reproduce the one-imidazole geometry and energy accurately, we hypothesize that they will reasonably predict the binding free energy of higher-level coordination numbers. Hence, we did not extend the analysis of these ions up to six imidazole complexes. Overall, the results shed light on metal-protein interactions by emphasizing the importance of ligand-ligand interaction and metal-π-stacking within metalloproteins.

A coordinative modular assembly-constructed self-reinforced nano-therapeutic agent for effective antitumor-immune activation

Immunosuppressive microenvironment and poor immunogenicity are two stumbling blocks in anti-tumor immune activation. Tumor associated macrophages (TAMs) play crucial roles in immunosuppressive microenvironment, while immunogenic cell death (ICD) is a typical strategy to boost immunogenicity. Herein, we developed a coordinative modular assembly-based self-reinforced nanoparticle, (CaO2/TA)-(Fe3+/BSA) which integrated CaO2, Fe3+-tannic acid coordinated networks and albumin under the instruction of molecular dynamics simulation. (CaO2/TA)-(Fe3+/BSA) could significantly enhance Fenton reaction through Fe3+ self-reduction and H2O2 self-sufficiency, and simultaneously increased intracellular accumulation of Ca2+. The self-augmented Fenton reaction with sufficient reactive oxygen species effectively repolarized TAMs and elicited ICD with Ca2+ overload. Besides, (CaO2/TA)-(Fe3+/BSA) was confirmed to self-reinforce deep tumor drug delivery by "treatment-delivery" positive feedback based on gp60-mediated transcytosis and M2-like macrophages repolarization-mediated perfusion promotion. Resultantly, (CaO2/TA)-(Fe3+/BSA) effectively alleviated immunosuppression, provoked local and systemic immune response and potentiated anti-PD-1 antibody therapy. Our strategy highlights a facile and controllable approach to construct penetrated effective antitumor nano-immunotherapeutic agent.

Binding of Phosphate Species to Ca2+ and Mg2+ in Aqueous Solution

Phosphate derivatives and their interaction with metal cations are involved in many important biological phenomena, so an accurate characterization of the phosphate-metal interaction is necessary to properly understand the role of phosphate-metal contacts in mediating biological function. Herein, we improved the standard 12-6 Lennard-Jones (LJ) potential via the usage of the 12-6-4 LJ model, which incorporates ion-induced dipole interactions. Via parameter scanning, we fine-tuned the 12-6-4 LJ polarizability values to obtain accurate absolute binding free energies for the phosphate anions H2PO4-, HPO42-, PO43- coordinating with Ca2+ and Mg2+. First, we modified the phosphate 12-6-4 LJ parameters to reproduce the solvation free energies of the series of phosphate anions using the thermodynamic integration (TI) method. Then, using the potential mean force (PMF) method, the polarizability of the metal-phosphate interaction was obtained. We show that the free energy profiles of phosphate ions coordinated to Ca2+ and Mg2+ generally show similar trends at longer metal-phosphate distances, while the absolute binding energy values increased with deprotonation. The resulting parameters demonstrate the flexibility of the 12-6-4 LJ-type nonbonded model and its usefulness in accurately describing cation-anion interactions.

Evidence of Sharp Transitions between Octahedral and Capped Trigonal Prism States of the Solvation Shell of the Aqueous Fe3+ Ion

The structure of the solvation shell of the aqueous Fe3+ ion has been a subject of controversy due to discrepancies between experiments and different levels of theory. We address this issue by performing simulations for a wide range of ion concentrations, using various potential energy functions, supplemented by density functional theory calculations of selected configurations. The solvation shell undergoes abrupt transitions between two states: a hexacoordinated octahedral (OH) state and a capped trigonal prism (CTP) state with 7-fold coordination. The lifetime of these states is dependent on concentration. In dilute FeCl3 solutions, the lifetimes of both are similar (≈1 ns). However, the lifetime of the OH state increases with ion concentration, while that of the CTP state decreases slightly. When a uniform negative background charge is used instead of explicit counterions, the lifetime of the OH state is greatly overestimated. These findings underscore the need for further experimental measurements and higher-level simulations.

Investigation of the insecticidal potential of curcumin derivatives that target the Helicoverpa armigera sterol carrier protein-2

Cotton bollworm (Helicoverpa armigera) is a highly polyphagous, widely prevalent, and persistent Old World insect pest that affects numerous important crops that are directly consumed by people, including tomato, cotton, pigeon pea, chickpea, rice, sorghum, and cowpea. Insects do not synthesize steroids but obtain them from their diet. Inhibition of dietary uptake of steroids by insects is a potentially effective insecticidal mechanism that should not be toxic to humans and other mammals, who synthesize their steroids. Ten curcumin derivatives were tested against H. armigera sterol carrier protein-2 (HaSCP-2) for their potential as insecticidal agents. Curcumin derivatives were initially docked at the binding site of HaSCP-2 to determine their binding affinities and plausible binding modes. The binding modes predominantly show hydrophobic interactions of derivatives with Phe53, Phe110, and Phe89 as core interacting residues in the active site. Validation of in silico results was carried out by performing a fluorescence binding and displacement assay to determine the binding affinities of curcumin derivatives. Among a collection of curcumin derivatives tested, Cur10 showed the lowest IC50 value of 9.64 μM, while Cur07 was 19.86 μM, and Cur06 was 20.79 μM. There was a significant negative correlation between the ability of the curcumin derivatives tested to displace the fluorescent probe from the sterol binding site of HaSCP-2 and to inhibit Sf9 insect cell growth in culture, which is consistent with the curcumin derivatives acting by the novel mechanism of blocking sterol uptake. Then molecular dynamics simulation studies validated the predicted binding modes and the interactions of curcumin derivatives with HaSCP-2 protein. In conclusion, these studies support the potential use of curcumin derivatives as insecticidal agents.

Multifunctional metal-coordinated Co-assembled carrier-free nanoplatform based on dual-drugs for ferroptosis-mediated cocktail therapy of hepatocellular carcinoma growth and metastasis

The heterogeneity of hepatocellular carcinoma (HCC) and the complexity of the tumor microenvironment (TME) pose challenges to efficient drug delivery and the antitumor efficacy of combined or synergistic therapies. Herein, a metal-coordinated carrier-free nanodrug (named as USFe3+ LA NPs) was developed for ferroptosis-mediated multimodal synergistic anti-HCC. Natural product ursolic acid (UA) was incorporated to enhance the sensitivity of tumor cells to sorafenib (SRF). Surface decoration of cell penetration peptide and epithelial cell adhesion molecule aptamer facilitated the uptake of USFe3+ LA NPs by HepG2 cells. Meanwhile, Fe3+ ions could react with intracellular hydrogen peroxide, generating toxic hydroxyl radical (·OH) for chemodynamical therapy (CDT) and amplified ferroptosis by cystine/glutamate antiporter system (System Xc-), which promoted the consumption of glutathione (GSH) and inhibited the expression of glutathione peroxidase 4 (GPX4). Notably, these all-in-one nanodrugs could inhibit tumor metastasis and induced immunogenic cell death (ICD). Last but not least, the nanodrugs demonstrated favorable biocompatibility, augmenting the immune response against the programmed death-ligand 1 (PD-L1) by increasing cytotoxic T cell infiltration. In vivo studies revealed significant suppression of tumor growth and distant metastasis. Overall, our work introduced a novel strategy for applications of metal-coordinated co-assembled carrier-free nano-delivery system in HCC combination therapy, especially in the realms of cancer metastasis prevention and immunotherapy.

The protein phosphorylation landscape in photosystem I of the desert algae Chlorella sp

The phosphorylation of photosystem II (PSII) and its antenna (LHCII) proteins has been studied, and its involvement in state transitions and PSII repair is known. Yet, little is known about the phosphorylation of photosystem I (PSI) and its antenna (LHCI) proteins. Here, we applied proteomics analysis to generate a map of the phosphorylation sites of the PSI-LHCI proteins in Chlorella ohadii cells that were grown under low or extreme high-light intensities (LL and HL). Furthermore, we analyzed the content of oxidized tryptophans and PSI-LHCI protein degradation products in these cells, to estimate the light-induced damage to PSI-LHCI. Our work revealed the phosphorylation of 17 of 22 PSI-LHCI subunits. The analyses detected the extensive phosphorylation of the LHCI subunits Lhca6 and Lhca7, which is modulated by growth light intensity. Other PSI-LHCI subunits were phosphorylated to a lesser extent, including PsaE, where molecular dynamic simulation proposed that a phosphoserine stabilizes ferredoxin binding. Additionally, we show that HL-grown cells accumulate less oxidative damage and degradation products of PSI-LHCI proteins, compared with LL-grown cells. The significant phosphorylation of Lhca6 and Lhca7 at the interface with other LHCI subunits suggests a physiological role during photosynthesis, possibly by altering light-harvesting characteristics and binding of other subunits.

Lanthanide transport in angstrom-scale MoS2-based two-dimensional channels

Rare earth elements (REEs), critical to modern industry, are difficult to separate and purify, given their similar physicochemical properties originating from the lanthanide contraction. Here, we systematically study the transport of lanthanide ions (Ln3+) in artificially confined angstrom-scale two-dimensional channels using MoS2-based building blocks in an aqueous environment. The results show that the uptake and permeability of Ln3+ assume a well-defined volcano shape peaked at Sm3+. This transport behavior is rooted from the tradeoff between the barrier for dehydration and the strength of interactions of lanthanide ions in the confinement channels, reminiscent of the Sabatier principle. Molecular dynamics simulations reveal that Sm3+, with moderate hydration free energy and intermediate affinity for channel interaction, exhibit the smallest dehydration degree, consequently resulting in the highest permeability. Our work not only highlights the distinct mass transport properties under extreme confinement but also demonstrates the potential of dialing confinement dimension and chemistry for greener REEs separation.

Development of a Fluctuating Charge Model for Zinc-Containing Metalloproteins

Metalloproteins widely exist in biology and play important roles in various processes. To accurately simulate metalloprotein systems, modeling polarization and charge transfer effects is vital. The fluctuating charge (FQ) model can efficiently generate atomic charges and simulate the charge transfer effect; it has been developed for a wide range of applications, but few models have been specifically tailored for metalloproteins. In this study, we present a fluctuating charge model specifically for zinc-containing metalloproteins based on the extended charge equilibration (EQeq) scheme. Our model was parametrized to reproduce CM5 charges instead of RESP/CHELPG charges because the former is less dependent on the conformation or basis set, does not suffer from unphysical charges for buried atoms, and is still being able to well reproduce the molecular dipoles. During our study, we found that adding the Pauling-bond-order-like term (referred to as the "+C term" in a previous study) between the zinc ion and ligating atoms significantly improves the model's performance. Although our model was trained for four-coordinated zinc sites only, our results indicated it can well describe the atomic charges in diverse zinc sites. Morever, our model was used to generate partial charges for the metal sites in three different zinc-containing metalloproteins (with four-, five-, and six-coordinated metal sites, respectively). These charges exhibited performance comparable to that of the RESP charges in molecular dynamics (MD) simulations. Additional tests indicated our model could also well reproduce the CM5 charges when geometric changes were involved. Those results indicate that our model can efficiently calculate the atomic charges for metal sites and well simulate the charge transfer effect, which marks an important step toward developing versatile polarizable force fields for metalloproteins.

Iron uptake pathway of Escherichia coli as an entry route for peptide nucleic acids conjugated with a siderophore mimic

Bacteria secrete various iron-chelators (siderophores), which scavenge Fe3+ from the environment, bind it with high affinity, and retrieve it inside the cell. After the Fe3+ uptake, bacteria extract the soluble iron(II) from the siderophore. Ferric siderophores are transported inside the cell via the TonB-dependent receptor system. Importantly, siderophore uptake paths have been also used by sideromycins, natural antibiotics. Our goal is to hijack the transport system for hydroxamate-type siderophores to deliver peptide nucleic acid oligomers into Escherichia coli cells. As siderophore mimics we designed and synthesized linear and cyclic Nδ-acetyl-Nδ-hydroxy-l-ornithine based peptides. Using circular dichroism spectroscopy, we found that iron(III) is coordinated by the linear trimer with hydroxamate groups but not by the cyclic peptide. The internal flexibility of the linear siderophore oxygen atoms and their interactions with Fe3+ were confirmed by all-atom molecular dynamics simulations. Using flow cytometry we found that the designed hydroxamate trimer transports PNA oligomers inside the E. coli cells. Growth recovery assays on various E. coli mutants suggest the pathway of this transport through the FhuE outer-membrane receptor, which is responsible for the uptake of the natural iron chelator, ferric-coprogen. This pathway also involves the FhuD periplasmic binding protein. Docking of the siderophores to the FhuE and FhuD receptor structures showed that binding of the hydroxamate trimer is energetically favorable corroborating the experimentally suggested uptake path. Therefore, this siderophore mimic, as well as its conjugate with PNA, is most probably internalized through the hydroxamate pathway.

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.

Structure and properties variations of regenerated cellulose fibers induced by metal ion impurity

In the ionic liquids (ILs) method for processing regenerated cellulose fiber (RCF), which is a high-performance ecologically benign product, metal ion impurities (such as Fe3+ and Cu2+) of cellulose might inevitably remain in the recycled ILs and coagulation bath. The presence of metal ions might negatively impact the properties of the manufactured RCFs and obstruct their applications, which are urgent to be made clear. For this research, the solvent for dissolving wood pulp cellulose (WPC) was 1-ethyl-3-methylimidazolium diethyl phosphate ([Emim]DEP) with various metal ion concentrations. The effect of metal ions in IL on the dissolution of cellulose was investigated by Molecular Dynamics simulations. Rheological analysis and degree of polymerization (DP) analysis were adopted to evaluate the influence on fiber spinnability of different spinning solution metal ion concentrations and various dissolving times. Further, the morphology and mechanical performances of the RCFs variation regulation were also thoroughly researched. The findings showed that the presence of metal ions in the spinning solution affected the DP, crystallinity, and orientation factor of RCFs, which will influence their stress more sensitively than the strain. These findings can serve as a practical guide for the commercial manufacture of emerging fiber.

Development of AMOEBA Polarizable Force Field for Rare-Earth La3+ Interaction with Bioinspired Ligands

Rare-earth metals (REMs) are crucial for many important industries, such as power generation and storage, in addition to cancer treatment and medical imaging. One promising new REM refinement approach involves mimicking the highly selective and efficient binding of REMs observed in relatively recently discovered proteins. However, realizing any such bioinspired approach requires an understanding of the biological recognition mechanisms. Here, we developed a new classical polarizable force field based on the AMOEBA framework for modeling a lanthanum ion (La3+) interacting with water, acetate, and acetamide, which have been found to coordinate the ion in proteins. The parameters were derived by comparing to high-level ab initio quantum mechanical (QM) calculations that include relativistic effects. The AMOEBA model, with advanced atomic multipoles and electronic polarization, is successful in capturing both the QM distance-dependent La3+-ligand interaction energies and experimental hydration free energy. A new scheme for pairwise polarization damping (POLPAIR) was developed to describe the polarization energy in La3+ interactions with both charged and neutral ligands. Simulations of La3+ in water showed water coordination numbers and ion-water distances consistent with previous experimental and theoretical findings. Water residence time analysis revealed both fast and slow kinetics in water exchange around the ion. This new model will allow investigation of fully solvated lanthanum ion-protein systems using GPU-accelerated dynamics simulations to gain insights on binding selectivity, which may be applied to the design of synthetic analogues.

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.

Simulating core electron binding energies of halogenated species adsorbed on ice surfaces and in solution via relativistic quantum embedding calculations

In this work, we investigate the effects of the environment on the X-ray photoelectron spectra of hydrogen chloride and chloride ions adsorbed on ice surfaces, as well as of chloride ions in water droplets. In our approach, we combine a density functional theory (DFT) description of the ice surface with that of halogen species using the recently developed relativistic core-valence separation equation of motion coupled cluster (CVS-EOM-IP-CCSD) via the frozen density embedding formalism (FDE), to determine the K and L_1,2,3 edges of chlorine. Our calculations, which incorporate temperature effects through snapshots from classical molecular dynamics simulations, are shown to reproduce experimental trends in the change of the core binding energies of Cl^- upon moving from a liquid (water droplets) to an interfacial (ice quasi-liquid layer) environment. Our simulations yield water valence band binding energies in good agreement with experiment, which vary little between the droplets and the ice surface. For halide core binding energies there is an overall trend for overestimating experimental values, though good agreement between theory and experiment is found for Cl^- in water droplets and on ice. For HCl on the other hand there are significant discrepancies between experimental and calculated core binding energies when we consider structural models that maintain the H-Cl bond more or less intact. An analysis of models that allow for pre-dissociated and dissociated structures suggests that experimentally observed chemical shifts in binding energies between Cl^- and HCl would require that H⁺ (in the form of H₃O⁺) and Cl^- are separated by roughly 4-6 Å.

Force Field Parameterization of Actinyl Molecular Cations Using the 12-6-4 Model

In this work, a set of 12-6-4 force fields (FFs) parameters were developed for the actinyl molecular cations, AnO2n+ (n = 1, 2), from uranium to plutonium for classical molecular dynamics (MD) for four water models: TIP3P, SPC/E, OPC3, and TIP4Pew. Such a non-bonded potential model taking into account the induced dipole between the metallic center and the surrounding molecules has shown better performances for various cations than the classic 12-6 non-bonded potentials. The parametrization method proposed elsewhere for metallic cations has been extended to these molecular cations. In contrast to the actinyl 12-6 FFs from the literature, the new models reproduce correctly both solvation and thermodynamic properties, thanks to the inclusion of the induced dipole term (C4). The transferability of such force fields was assessed by performing MD simulations of carbonato actinyl species, which are highly implicated in actinide migration or actinide extraction from seawater. A highly satisfying agreement was found when comparing the EXAFS signals computed from our MD simulation to the experimental ones. The set of FFs developed here opens new possibilities for the study of actinide chemistry.

Systematic Evaluation of Ion Diffusion and Water Exchange

As a fundamental property of all fluids, diffusion plays myriad roles in both science and our daily lives. Diffusive properties of many liquids including water have been extensively studied both experimentally and theoretically, while for transition metal ions, there exist significant experimental data that have not been extensively studied theoretically. Hence, high-confidence predictions for challenging systems like radioactive ions that are biohazardous cannot be reliably generated. In this work, a workflow named ISAIAH (Ion Simulation using AMBER for dIffusion Action when Hydrated) was designed to accurately simulate the diffusion coefficients of 15 monoatomic ions with charges varying from -1 to +3 in four water models. As the results indicate, good agreement with experimental values was achieved, leading us to select ²³⁹Pu⁴⁺ (for which no experimental data are available) as a candidate ion to make a theoretical prediction of its diffusion coefficient in water. Among all the force field parameter sets, the ones parametrized using an augmented 12-6-4 Lennard-Jones (LJ) potential showed lower average unsigned errors (AUE) for ions of various radii and electron configurations relative to some 12-6 LJ parameters. This observation agrees well with the fact that diffusion is affected by both the hydration free energy (HFE) and the ion-oxygen distance (IOD) between solute and solvent molecules, both of which are handled well by the 12-6-4 model.

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.

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.

Bridging the 12-6-4 Model and the Fluctuating Charge Model

Metal ions play important roles in various biological systems. Molecular dynamics (MD) using classical force field has become a popular research tool to study biological systems at the atomic level. However, meaningful MD simulations require reliable models and parameters. Previously we showed that the 12-6 Lennard-Jones nonbonded model for ions could not reproduce the experimental hydration free energy (HFE) and ion-oxygen distance (IOD) values simultaneously when ion has a charge of +2 or higher. We discussed that this deficiency arises from the overlook of the ion-induced dipole interaction in the 12-6 model, and this term is proportional to 1/r⁴ based on theory. Hence, we developed the 12-6-4 model and showed it could solve this deficiency in a physically meaningful way. However, our previous research also found that the 12-6-4 model overestimated the coordination numbers (CNs) for some highly charged metal ions. And we attributed this artifact to that the current 12-6-4 scheme lacks a correction for the interactions among the first solvation shell water molecules. In the present study, we considered the ion-included dipole interaction by using the 12-6 model with adjusting the atomic charges of the first solvation shell water molecules. This strategy not only considers the ion-induced dipole interaction between ion and the first solvation shell water molecules but also well accounts for the increased repulsion among these water molecules compared to the bulk water molecules. We showed this strategy could well reproduce the experimental HFE and IOD values for Mg²⁺, Zn²⁺, Al³⁺, Fe³⁺, and In³⁺ and solve the CN overestimation issue of the 12-6-4 model for Fe³⁺ and In³⁺. Moreover, our simulation results showed good agreement with previous ab initio MD simulations. In addition, we derived the physical relationship between the C₄ parameter and induced dipole moment, which agreed well with our simulation results. Finally, we discussed the implications of the present work for simulating metalloproteins. Due to the fluctuating charge model uses a similar concept to the 12-6 model with adjusting atomic charges, we believe the present study builds a bridge between the 12-6-4 model and the fluctuating charge model.

Advanced Electrostatic Model for Monovalent Ions Based on Ab Initio Energy Decomposition

Ions play important roles in the structures and functions of biomolecules. In biomolecular simulations, ions either directly interact with biomolecules or provide an ionic environment that influences electrostatic interactions of solutes. The AMOEBA+ water model has demonstrated significant advancement of the classical force field for describing molecular interactions due to its improvements on the functional forms to account for essential physics. This work expands the applicability of the AMOEBA+ model toward alkali metal (Li, Na, K, Rb, and Cs) and halogen (F, Cl, Br, and I) ions. Various quantum chemical data on ion-ion and ion-water interactions, experimental ion hydration free energies, and lattice energies of salt crystals are used in the parametrization. The final parameters are verified with other properties outside of the parametrization data, including lattice energies of additional salt crystals and ionic activity coefficients in solution. The new model captures a wide range of ion properties from the gas phase to solution phase and crystals. More importantly, AMOEBA+ provides energy components that are consistent with ab initio energy decomposition. Thus, we expect AMOEBA+ to be more general, transferable, and valuable for the interpretation of intermolecular forces in efficient classical simulations.

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.