Rational Design of Methodology-Independent Metal Parameters Using a Nonbonded Dummy Model

Computational Modeling Study of the Molecular Basis of dNTP Selectivity in Human Terminal Deoxynucleotidyltransferase

Human terminal deoxynucleotidyl transferase (TdT) can catalyze template-independent DNA synthesis during the V(D)J recombination and DNA repair through nonhomologous end joining. The capacity for template-independent random addition of nucleotides to single-stranded DNA makes this polymerase useful in various molecular biological applications involving sequential stepwise synthesis of oligonucleotides using modified dNTP. Nonetheless, a serious limitation to the applications of this enzyme is strong selectivity of human TdT toward dNTPs in the order dGTP > dTTP ≈ dATP > dCTP. This study involved molecular dynamics to simulate a potential impact of amino acid substitutions on the enzyme's selectivity toward dNTPs. It was found that the formation of stable hydrogen bonds between a nitrogenous base and amino acid residues at positions 395 and 456 is crucial for the preferences for dNTPs. A set of single-substitution and double-substitution mutants at these positions was analyzed by molecular dynamics simulations. The data revealed two TdT mutants-containing either substitution D395N or substitutions D395N+E456N-that possess substantially equalized selectivity toward various dNTPs as compared to the wild-type enzyme. These results will enable rational design of TdT-like enzymes with equalized dNTP selectivity for biotechnological applications.

Computational design of carboxylase for the synthesis of 4-hydroxyisophthalic acid from p-hydroxybenzoic acid by fixing CO2

Carbon dioxide (CO2) emissions constitute the primary contribution to global climate change. Synthetic CO2 fixation represents an exceptionally appealing and sustainable method for carbon neutralization. Unlike the limitations of chemical catalysis, biological CO2 fixation displays high selectivity and the ability to operate under mild conditions. The superfamily of amidohydrolases has demonstrated the ability to synthesize a range of aromatic monocarboxylic acids. However, there is a scarcity of reported carboxylases capable of synthesizing aromatic dicarboxylic acids. Among these, 4-hydroxyisophthalic acid holds significant potential for applications across various fields, yet no enzyme has been reported for its synthesis. In this study, we developed for the first time that exhibits starting activity in fixing CO2 to synthesize 4-hydroxyisophthalic acid. Furthermore, we have devised a computational strategy that effectively enhances the catalytic activity of this enzyme. A focused library comprising only 13 variants was generated. Experimental validation confirmed a threefold improvement in the carboxylation activity of the optimal variant (L47M). The computational enzyme design strategy proposed in this paper demonstrates broad applicability in developing carboxylases for synthesizing other aromatic dicarboxylic acids. This lays the groundwork for leveraging biocatalysis in industrial synthesis for CO2 fixation.

SNP-Associated Substitutions of Amino Acid Residues in the dNTP Selection Subdomain Decrease Polβ Polymerase Activity

In the cell, DNA polymerase β (Polβ) is involved in many processes aimed at maintaining genome stability and is considered the main repair DNA polymerase participating in base excision repair (BER). Polβ can fill DNA gaps formed by other DNA repair enzymes. Single-nucleotide polymorphisms (SNPs) in the POLB gene can affect the enzymatic properties of the resulting protein, owing to possible amino acid substitutions. For many SNP-associated Polβ variants, an association with cancer, owing to changes in polymerase activity and fidelity, has been shown. In this work, kinetic analyses and molecular dynamics simulations were used to examine the activity of naturally occurring polymorphic variants G274R, G290C, and R333W. Previously, the amino acid substitutions at these positions have been found in various types of tumors, implying a specific role of Gly-274, Gly-290, and Arg-333 in Polβ functioning. All three polymorphic variants had reduced polymerase activity. Two substitutions-G274R and R333W-led to the almost complete disappearance of gap-filling and primer elongation activities, a decrease in the deoxynucleotide triphosphate-binding ability, and a lower polymerization constant, due to alterations of local contacts near the replaced amino acid residues. Thus, variants G274R, G290C, and R333W may be implicated in an elevated level of unrepaired DNA damage.

The Impact of SNP-Induced Amino Acid Substitutions L19P and G66R in the dRP-Lyase Domain of Human DNA Polymerase β on Enzyme Activities

Base excision repair (BER), which involves the sequential activity of DNA glycosylases, apurinic/apyrimidinic endonucleases, DNA polymerases, and DNA ligases, is one of the enzymatic systems that preserve the integrity of the genome. Normal BER is effective, but due to single-nucleotide polymorphisms (SNPs), the enzymes themselves-whose main function is to identify and eliminate damaged bases-can undergo amino acid changes. One of the enzymes in BER is DNA polymerase β (Polβ), whose function is to fill gaps in DNA. SNPs can significantly affect the catalytic activity of an enzyme by causing an amino acid substitution. In this work, pre-steady-state kinetic analyses and molecular dynamics simulations were used to examine the activity of naturally occurring variants of Polβ that have the substitutions L19P and G66R in the dRP-lyase domain. Despite the substantial distance between the dRP-lyase domain and the nucleotidyltransferase active site, it was found that the capacity to form a complex with DNA and with an incoming dNTP is significantly altered by these substitutions. Therefore, the lower activity of the tested polymorphic variants may be associated with a greater number of unrepaired DNA lesions.

A transferable double exponential potential for condensed phase simulations of small molecules

The Lennard-Jones potential is the most widely-used function for the description of non-bonded interactions in transferable force fields for the condensed phase. This is not because it has an optimal functional form, but rather it is a legacy resulting from when computational expense was a major consideration and this potential was particularly convenient numerically. At present, it persists because the effort that would be required to re-write molecular modelling software and train new force fields has, until now, been prohibitive. Here, we present Smirnoff-plugins as a flexible framework to extend the Open Force Field software stack to allow custom force field functional forms. We deploy Smirnoff-plugins with the automated Open Force Field infrastructure to train a transferable, small molecule force field based on the recently-proposed double exponential functional form, on over 1000 experimental condensed phase properties. Extensive testing of the resulting force field shows improvements in transfer free energies, with acceptable conformational energetics, run times and convergence properties compared to state-of-the-art Lennard-Jones based force fields.

Individual Contributions of Amido Acid Residues Tyr122, Ile168, and Asp173 to the Activity and Substrate Specificity of Human DNA Dioxygenase ABH2

Human Fe(II)/α-ketoglutarate-dependent dioxygenase ABH2 plays a crucial role in the direct reversal repair of nonbulky alkyl lesions in DNA nucleobases, e.g., N¹-methyladenine (m¹A), N³-methylcytosine (m³C), and some etheno derivatives. Moreover, ABH2 is capable of a less efficient oxidation of an epigenetic DNA mark called 5-methylcytosine (m⁵C), which typically is a specific target of DNA dioxygenases from the TET family. In this study, to elucidate the mechanism of the substrate specificity of ABH2, we investigated the role of several active-site amino acid residues. Functional mapping of the lesion-binding pocket was performed through the analysis of the functions of Tyr122, Ile168, and Asp173 in the damaged base recognition mechanism. Interactions of wild-type ABH2, or its mutants Y122A, I168A, or D173A, with damaged DNA containing the methylated base m¹A or m³C or the epigenetic marker m⁵C were analyzed by molecular dynamics simulations and kinetic assays. Comparative analysis of the enzymes revealed an effect of the substitutions on DNA binding and on catalytic activity. Obtained data clearly demonstrate the effect of the tested amino acid residues on the catalytic activity of the enzymes rather than the DNA-binding ability. Taken together, these data shed light on the molecular and kinetic consequences of the substitution of active-site residues for the mechanism of the substrate recognition.

The Activity of Natural Polymorphic Variants of Human DNA Polymerase β Having an Amino Acid Substitution in the Transferase Domain

To maintain the integrity of the genome, there is a set of enzymatic systems, one of which is base excision repair (BER), which includes sequential action of DNA glycosylases, apurinic/apyrimidinic endonucleases, DNA polymerases, and DNA ligases. Normally, BER works efficiently, but the enzymes themselves (whose primary function is the recognition and removal of damaged bases) are subject to amino acid substitutions owing to natural single-nucleotide polymorphisms (SNPs). One of the enzymes in BER is DNA polymerase β (Polβ), whose function is to fill gaps in DNA with complementary dNMPs. It is known that many SNPs can cause an amino acid substitution in this enzyme and a significant decrease in the enzymatic activity. In this study, the activity of four natural variants of Polβ, containing substitution E154A, G189D, M236T, or R254I in the transferase domain, was analyzed using molecular dynamics simulations and pre-steady-state kinetic analyses. It was shown that all tested substitutions lead to a significant reduction in the ability to form a complex with DNA and with incoming dNTP. The G189D substitution also diminished Polβ catalytic activity. Thus, a decrease in the activity of studied mutant forms may be associated with an increased risk of damage to the genome.

Benchmarking biomolecular force field-based Zn2+ for mono- and bimetallic ligand binding sites

Zn²⁺ is one of the most versatile biologically available metal ions, but accurate modeling of Zn²⁺ -containing metalloproteins at the biomolecular force field level can be challenging. Since most Zn²⁺ models are parameterized in bulk solvent, in-depth knowledge about their performance in a protein environment is limited. Thus, we systematically investigate here the behavior of non-polarizable Zn²⁺ models for their ability to reproduce experimentally determined metal coordination and ligand binding in metalloproteins. The benchmarking is performed in challenging environments, including mono- (carbonic anhydrase II) and bimetallic (metallo-β-lactamase VIM-2) ligand binding sites. We identify key differences in the performance between the Zn²⁺ models with regard to the preferred ligating atoms (charged/non-charged), attraction of water molecules, and the preferred coordination geometry. Based on these results, we suggest suitable simulation conditions for varying Zn²⁺ site geometries that could guide the further development of biomolecular Zn²⁺ models.

Comparison of Empirical Zn2+ Models in Protein–DNA Complexes

Zinc ions are the second most abundant ions found in humans. Their role in proteins can be merely structural but also catalytic, owing to their transition metal character. Modelling their geometric–coordination versatility by empirical force fields is, thus, a challenging task. In this work, we evaluated three popular models, specifically designed to represent zinc ions with regard to their capability of preserving structural integrity. To this end, we performed molecular dynamics simulations of two zinc-containing protein–DNA complexes, which differed in their zinc coordination, i.e., four cysteines or two cysteines and two histidines. The most flexible non-bonded 12-6-4 Lennard–Jones-type model shows a preference for six-fold coordination of the Zn2+-ions in contradiction to the crystal structure. The cationic dummy atom model favours tetrahedral geometry, whereas the bonded extended zinc AMBER force field model, by construction, best preserves the initial geometry of a regular or slightly distorted tetrahedron. Our data renders the extended zinc AMBER force field the best model for structural zinc ions in a given geometry. In more complicated cases, though, more flexible models may be advantageous.

Human Polβ Natural Polymorphic Variants G118V and R149I Affects Substate Binding and Catalysis

DNA polymerase β (Polβ) expression is essential for the cell's response to DNA damage that occurs during natural cellular processes. Polβ is considered the main reparative DNA polymerase, whose role is to fill the DNA gaps arising in the base excision repair pathway. Mutations in Polβ can lead to cancer, neurodegenerative diseases, or premature aging. Many single-nucleotide polymorphisms have been identified in the POLB gene, but the consequences of these polymorphisms are not always clear. It is known that some polymorphic variants in the Polβ sequence reduce the efficiency of DNA repair, thereby raising the frequency of mutations in the genome. In the current work, we studied two polymorphic variants (G118V and R149I separately) of human Polβ that affect its DNA-binding region. It was found that each amino acid substitution alters Polβ's affinity for gapped DNA. Each polymorphic variant also weakens its binding affinity for dATP. The G118V variant was found to greatly affect Polβ's ability to fill gapped DNA and slowed the catalytic rate as compared to the wild-type enzyme. Thus, these polymorphic variants seem to decrease the ability of Polβ to maintain base excision repair efficiency.

Insight into the mechanism of DNA synthesis by human terminal deoxynucleotidyltransferase

Terminal deoxynucleotidyltransferase (TdT) is a member of the DNA polymerase X family that is responsible for random addition of nucleotides to single-stranded DNA. We present investigation into the role of metal ions and specific interactions of dNTP with active-site amino acid residues in the mechanisms underlying the recognition of nucleoside triphosphates by human TdT under pre-steady-state conditions. In the elongation mode, the ratios of translocation and dissociation rate constants, as well as the catalytic rate constant were dependent on the nature of the nucleobase. Preferences of TdT in dNTP incorporation were researched by molecular dynamics simulations of complexes of TdT with a primer and dNTP or with the elongated primer. Purine nucleotides lost the "summarised" H-bonding network after the attachment of the nucleotide to the primer, whereas pyrimidine nucleotides increased the number and relative lifetime of H-bonds in the post-catalytic complex. The effect of divalent metal ions on the primer elongation revealed that Me²⁺ cofactor can significantly change parameters of the primer elongation by strongly affecting the rate of nucleotide attachment and the polymerisation mode.

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.

Unlocking the computational design of metal-organic cages

Metal-organic cages are macrocyclic structures that can possess an intrinsic void that can hold molecules for encapsulation, adsorption, sensing, and catalysis applications. As metal-organic cages may be comprised from nearly any combination of organic and metal-containing components, cages can form with diverse shapes and sizes, allowing for tuning toward targeted properties. Therefore, their near-infinite design space is almost impossible to explore through experimentation alone and computational design can play a crucial role in exploring new systems. Although high-throughput computational design and screening workflows have long been known as powerful tools in drug and materials discovery, their application in exploring metal-organic cages is more recent. We show examples of structure prediction and host-guest/catalytic property evaluation of metal-organic cages. These examples are facilitated by advances in methods that handle metal-containing systems with improved accuracy and are the beginning of the development of automated cage design workflows. We finally outline a scope for how high-throughput computational methods can assist and drive experimental decisions as the field pushes toward functional and complex metal-organic cages. In particular, we highlight the importance of considering realistic, flexible systems.

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.

Ribosome Elongation Kinetics of Consecutively Charged Residues Are Coupled to Electrostatic Force

The speed of protein synthesis can dramatically change when consecutively charged residues are incorporated into an elongating nascent protein by the ribosome. The molecular origins of this class of allosteric coupling remain unknown. We demonstrate, using multiscale simulations, that positively charged residues generate large forces that move the P-site amino acid away from the A-site amino acid. Negatively charged residues generate forces of similar magnitude but move the A- and P-sites closer together. These conformational changes, respectively, increase and decrease the transition state barrier height to peptide bond formation, explaining how charged residues mechanochemically alter translation speed. This mechanochemical mechanism is consistent with in vivo ribosome profiling data exhibiting proportionality between translation speed and the number of charged residues, experimental data characterizing nascent chain conformations, and a previously published cryo-EM structure of a ribosome-nascent chain complex containing consecutive lysines. These results expand the role of mechanochemistry in translation and provide a framework for interpreting experimental results on translation speed.

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.

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.

Mechanical Forces Have a Range of Effects on the Rate of Ribosome Catalyzed Peptidyl Transfer Depending on Direction

Mechanical forces acting on the nascent chain residue located at the P-site of the ribosome can influence codon translation rates. Most observations to date involve force vectors aligned collinear with the long axis of the ribosome exit tunnel. What is poorly understood is how force applied in other directions will impact the rate of peptide bond formation catalyzed by the ribosome. Here, we utilize quantum mechanical/molecular mechanics simulations to estimate the changes in the activation free energy as a consequence of applying a constant force in various directions on the C-terminal residue at the P-site. Qualitatively consistent with the Bell model, we find this force can either accelerate, decelerate, or not alter the reaction rate depending on the force direction. A force in the average direction between the P-site 3' O-C ester bond that breaks and the peptide bond that forms accelerates the reaction. A force in the opposite direction slows down the reaction as it opposes these bonds breaking and forming, but surprisingly it does not do so to the maximum extent possible. In this case, there is a counterbalancing trend; the force in this direction brings the A-site amino nitrogen and the P-site tRNA A76 3' oxygen groups closer together, which promotes one of the proton shuttling steps of the reaction. We find the maximum force-induced slowdown occurs 37° off this axis. If force is applied in orthogonal directions to the reaction coordinates, there is no significant change in the reaction speed. These results indicate that there is a richer set of scenarios of force effects on translation speed that have yet to be experimentally explored and raise the possibility that cells could use these mechanochemical effects to modulate and regulate protein synthesis.

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.

Determination of van der Waals Parameters Using a Double Exponential Potential for Nonbonded Divalent Metal Cations in TIP3P Solvent

A double exponential (DE) functional form for Lennard-Jones (LJ) interactions, proposed in our previous study, has many advantages over LJ potentials including a natural softcore characteristic for the convenience of the pathway-based free-energy calculations, fast convergence, and flexibility in use. In this work, we put the first step on the application of the DE functional form by identifying a DE potential, coined DE-TIP3P, for molecular simulations using the TIP3P water model. The developed DE-TIP3 potential was better than LJ potential in reproducing the experimental water properties. Afterward, we developed the nonbonded models of 15 divalent metal ions, which frequently appear and play vital roles in biological systems, to be consistent with the DE-TIP3P potential and TIP3P water model. Our nonbonded models were as good as the complicated nonbonded dummy cationic models by Jiang et al. and the nonbonded 12-6-4 LJ models by Li and Merz in reproducing the experimental properties of those ions. Moreover, our nonbonded models achieved a better performance than the compromise (CM) LJ models and 12-6-4 LJ models, developed by Li and Merz, in reproducing the properties of MgCl₂ in aqueous solution.

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

Divalent metal ions play important roles in biological and materials systems. Molecular dynamics simulation is an efficient tool to investigate these systems at the microscopic level. Recently, four new water models (OPC3, OPC, TIP3P-FB, and TIP4P-FB) have been developed and better represent the physical properties of water than previous models. Metal ion parameters are dependent on the water model employed, making it necessary to develop metal ion parameters for select new water models. In the present work, we performed parameter scanning for the 12-6 Lennard-Jones nonbonded model of divalent metal ions in conjunction with the four new water models as well as four previous water models (TIP3P, SPC/E, TIP4P, and TIP4P-Ew). We found that these new three-point and four-point water models provide comparable or significantly improved performance for the simulation of divalent metal ions when compared to previous water models in the same category. Among all eight water models, the OPC3 water model yields the best performance for the simulation of divalent metal ions in the aqueous phase when using the 12-6 model. On the basis of the scanning results, we independently parametrized the 12-6 model for 24 divalent metal ions with each of the four new water models. As noted previously, the 12-6 model still fails to simultaneously reproduce the experimental hydration free energy (HFE) and ion-oxygen distance (IOD) values even with these new water models. To solve this problem, we parametrized the 12-6-4 model for the 16 divalent metal ions for which we have both experimental HFE and IOD values for each of the four new water models. The final parameters are able to reproduce both the experimental HFE and IOD values accurately. To validate the transferability of our parameters, we carried out benchmark calculations to predict the energies and geometries of ion-water clusters as well as the ion diffusivity coefficient of Mg²⁺. By comparison to quantum chemical calculations and experimental data, these results show that our parameters are well designed and have excellent transferability. The metal ion parameters for the 12-6 and 12-6-4 models reported herein can be employed in simulations of various biological and materials systems when using the OPC3, OPC, TIP3P-FB, or TIP4P-FB water model.

Comparing Nonbonded Metal Ion Models in the Divalent Cation Binding Protein PsaA

Divalent metal cations are essential for many biological processes; however, accurately modeling divalent metal ions has proved a significant challenge for molecular dynamics force fields. Here we show that the choice of ion model influences the observed dynamics in PsaA, a metal binding protein from Streptococcus pneumoniae. We conduct extensive unbiased simulations and free energy calculations of PsaA bound to its cognate ligand Mn²⁺ and inhibitory ligand Zn²⁺ using three nonbonded ion models: a 12-6 model, a 12-6-4 model, and a multisite model. The observed coordination geometries and metal binding dynamics are sensitive to the choice of ion model, with the most dramatic differences observed in free energy calculations of ion release. We show that the conformational ensemble of Mn-bound PsaA is more similar to the crystallographic metal bound open state. This work extends the current model of PsaA metal binding and provides a framework for the rationalization of experimentally determined metal binding behavior. Our findings support the use of the 12-6-4 ion model for further simulations of divalent cation binding proteins.

Evidence for an egg-box-like structure in iron(ii)–polygalacturonate hydrogels: a combined EXAFS and molecular dynamics simulation study

The coordination of Fe(ii) with polygalacturonic acid (polyGalA) in Fe(ii)–polyGalA hydrogels exhibits an octahedral geometry that follows the “egg-box model”.

Assessing the Performance of the Nonbonded Mg2+ Models in a Two-Metal-Dependent Ribonuclease

Magnesium ions (Mg²⁺), abundant in living cells, are essential for biomolecular structure, dynamics, and function. The biological importance of Mg²⁺ has motivated continuous development and improvement of various Mg²⁺ models for molecular dynamics (MD) simulations during the last decades. There are four types of nonbonded Mg²⁺ models: the point charge models based on a 12-6 or 12-6-4 type Lennard-Jones (LJ) potential, and the multisite models based on a 12-6 or 12-6-4 LJ potential. Here, we systematically assessed the performance of these four types of nonbonded Mg²⁺ models (21 models in total) in terms of maintaining a challenging intermediate state configuration captured in the structure of a prototypical two-metal-ion RNase H complex with an RNA/DNA hybrid. Our data demonstrate that the 12-6-4 multisite models, which account for charge-induced dipole interactions, perform the best in reproducing all the unique coordination modes in this intermediate state and maintaining the correct carboxylate denticity. Our benchmark work provides a useful guideline for MD simulations and structural refinement of Mg²⁺-containing biomolecular systems.

Single-Ion Thermodynamics from First Principles: Calculation of the Absolute Hydration Free Energy and Single-Electrode Potential of Aqueous Li+ Using ab Initio Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulations

A recently proposed thermodynamic integration (TI) approach formulated in the framework of quantum mechanical/molecular mechanical molecular dynamics (QM/MM MD) simulations is applied to study the structure, dynamics, and absolute intrinsic hydration free energy Δ_s G_M⁺,wat^◦ of the Li⁺ ion at a correlated ab initio level of theory. Based on the results, standard values (298.15 K, ideal gas at 1 bar, ideal solute at 1 molal) for the absolute intrinsic hydration free energy [Formula: see text] of the proton, the surface electric potential jump χ_wat^◦ upon entering bulk water, and the absolute single-electrode potential [Formula: see text] of the reference hydrogen electrode are calculated to be -1099.9 ± 4.2 kJ·mol^-1, 0.13 ± 0.08 V, and 4.28 ± 0.04 V, respectively, in excellent agreement with the standard values recommended by Hünenberger and Reif on the basis of an extensive evaluation of the available experimental data (-1100 ± 5 kJ·mol^-1, 0.13 ± 0.10 V, and 4.28 ± 0.13 V). The simulation results for Li⁺ are also compared to those for Na⁺ and K⁺ from a previous study in terms of relative hydration free energies ΔΔ_s G_M⁺,wat^◦ and relative electrode potentials [Formula: see text]. The calculated values are found to agree extremely well with the experimental differences in standard conventional hydration free energies ΔΔ_s G_M⁺,wat^• and redox potentials [Formula: see text]. The level of agreement between simulation and experiment, which is quantitative within error bars, underlines the substantial accuracy improvement achieved by applying a highly demanding QM/MM approach at the resolution-of-identity second-order Møller-Plesset perturbation (RIMP2) level over calculations relying on purely molecular mechanical or density functional theory (DFT) descriptions. A detailed analysis of the structural and dynamical properties of the Li⁺ hydrate indicates that a correct description of the solvation structure and dynamics is achieved as well at this level of theory. Consideration of the QM/MM potential-energy components also shows that the partitioning into QM and MM zones does not induce any significant energetic artifact for the system considered.

RNA Structural Dynamics As Captured by Molecular Simulations: A Comprehensive Overview

With both catalytic and genetic functions, ribonucleic acid (RNA) is perhaps the most pluripotent chemical species in molecular biology, and its functions are intimately linked to its structure and dynamics. Computer simulations, and in particular atomistic molecular dynamics (MD), allow structural dynamics of biomolecular systems to be investigated with unprecedented temporal and spatial resolution. We here provide a comprehensive overview of the fast-developing field of MD simulations of RNA molecules. We begin with an in-depth, evaluatory coverage of the most fundamental methodological challenges that set the basis for the future development of the field, in particular, the current developments and inherent physical limitations of the atomistic force fields and the recent advances in a broad spectrum of enhanced sampling methods. We also survey the closely related field of coarse-grained modeling of RNA systems. After dealing with the methodological aspects, we provide an exhaustive overview of the available RNA simulation literature, ranging from studies of the smallest RNA oligonucleotides to investigations of the entire ribosome. Our review encompasses tetranucleotides, tetraloops, a number of small RNA motifs, A-helix RNA, kissing-loop complexes, the TAR RNA element, the decoding center and other important regions of the ribosome, as well as assorted others systems. Extended sections are devoted to RNA-ion interactions, ribozymes, riboswitches, and protein/RNA complexes. Our overview is written for as broad of an audience as possible, aiming to provide a much-needed interdisciplinary bridge between computation and experiment, together with a perspective on the future of the field.

Extending the Nonbonded Cationic Dummy Model to Account for Ion-Induced Dipole Interactions

Modeling metalloproteins often requires classical molecular dynamics (MD) simulations in order to capture their relevant motions, which in turn necessitates reliable descriptions of the metal centers involved. One of the most successful approaches to date is provided by the "cationic dummy model", where the positive charge of the metal ion is transferred toward dummy particles that are bonded to the central metal ion in a predefined coordination geometry. While this approach allows for ligand exchange, and captures the correct electrostatics as demonstrated for different divalent metal ions, current dummy models neglect ion-induced dipole interactions. In the present work, we resolve this weakness by taking advantage of the recently introduced 12-6-4 type Lennard-Jones potential to include ion-induced dipole interactions. We revise our previous dummy model for Mg²⁺ and demonstrate that the resulting model can simultaneously reproduce the experimental solvation free energy and metal-ligand distances without the need for artificial restraints or bonds. As ion-induced dipole interactions become particularly important for highly charged metal ions, we develop dummy models for the biologically relevant ions Al³⁺, Fe³⁺, and Cr³⁺. Finally, the effectiveness of our new models is demonstrated in MD simulations of several diverse (and highly challenging to simulate) metalloproteins.

Free-Energy Calculations of Ionic Hydration Consistent with the Experimental Hydration Free Energy of the Proton

Computational free-energy correction strategies and the choice of experimental proton hydration free energy, ΔG_s^*(H⁺), are analyzed to investigate the apparent controversy in experimental thermodynamics of ionic hydration. Without corrections, the hydration free-energy (ΔG_hyd) calculations match experiments with ΔG_s^*(H⁺) = -1064 kJ/mol as reference. Using the Galvani surface potential the resulting (real) ΔG_hyd are consistent with ΔG_s^*(H⁺) = -1098 kJ/mol. When applying, in an ad hoc manner, the discrete solvent correction, ΔG_hyd matching the "consensus" ΔG_s^*(H⁺) of -1112 kJ/mol are obtained. This analysis rationalizes reports on ΔG_hyd calculations for ions using different experimental references. For neutral amino acid side chains ΔG_hyd are independent of the water model, whereas there are large differences in ΔG_hyd due to the water model for charged species, suggesting that long-range ordering of water around ions yields an important contribution to the ΔG_hyd. These differences are reduced significantly when applying consistent corrections, but to obtain the most accurate results it is recommended to use the water model belonging to the force field.

Development of a multisite model for Ni(II) ion in solution from thermodynamic and kinetic data

Force-field parameters are developed for a multisite model of Ni(II) ions to be used in molecular dynamics simulations combined to enhanced sampling methods. The performances of two charge-partitioning schemes are validated by taking into account structural, thermodynamic, and kinetic observables. One of the two models, featuring partial charges on the dummy atoms only, matches both Ni(II) free energy of solvation and water exchange rates. Such model is particularly suited to study complexation events at a fully dynamic description. © 2017 Wiley Periodicals, Inc.

Modeling Coordination-Directed Self-Assembly of M2L4 Nanocapsule Featuring Competitive Guest Encapsulation

Exploring the mechanism of self-assembly and guest encapsulation of nanocapsules is highly imperative for the design of sophisticated molecular containers and multistimuli-responsive functional materials. Here we present a molecular dynamics simulation protocol with implicit solvent and simulated annealing techniques to investigate the self-assembly and competitive guest (C₆₀ and C₇₀ fullerenes) encapsulation of a M₂L₄ nanocapsule that is self-assembled by the coordination of mercury cations and bent bidentate ligands. Stepwise formation of the nanocapsule and competitive fullerene encapsulation during dynamic structural changes in the self-assembly were detected successfully. Such processes were driven by coordination bonding and π-π stacking and obey the minimum total potential energy principle. Potential of mean force calculations for guest binding to the M₂L₄ nanocapsule explained the mechanism underlying the competitive encapsulations of C₆₀ and C₇₀. This work helps design new functional nanomaterials capable of guest encapsulation and release.

Metal Ion Modeling Using Classical Mechanics

Metal ions play significant roles in numerous fields including chemistry, geochemistry, biochemistry, and materials science. With computational tools increasingly becoming important in chemical research, methods have emerged to effectively face the challenge of modeling metal ions in the gas, aqueous, and solid phases. Herein, we review both quantum and classical modeling strategies for metal ion-containing systems that have been developed over the past few decades. This Review focuses on classical metal ion modeling based on unpolarized models (including the nonbonded, bonded, cationic dummy atom, and combined models), polarizable models (e.g., the fluctuating charge, Drude oscillator, and the induced dipole models), the angular overlap model, and valence bond-based models. Quantum mechanical studies of metal ion-containing systems at the semiempirical, ab initio, and density functional levels of theory are reviewed as well with a particular focus on how these methods inform classical modeling efforts. Finally, conclusions and future prospects and directions are offered that will further enhance the classical modeling of metal ion-containing systems.

Functional roles of magnesium binding to extracellular signal-regulated kinase 2 explored by molecular dynamics simulations and principal component analysis

Molecular dynamics (MD) simulations coupled with principal component (PC) analysis were carried out to study functional roles of Mg²⁺ binding to extracellular signal-regulated kinase 2 (ERK2). The results suggest that Mg²⁺ binding heavily decreases eigenvalue of the first principal component and totally inhibits motion strength of ERK2, which favors stabilization of ERK2 structure. Binding free energy predictions indicate that Mg²⁺ binding produces an important effect on binding ability of adenosine triphosphate (ATP) to ERK2 and strengthens the ATP binding. The calculations of residue-based free energy decomposition show that lack of Mg²⁺ weakens interactions between the hydrophobic rings of ATP and five residues I29, V37, A50, L105, and L154. Hydrogen bond analyses also prove that Mg²⁺ binding increases occupancies of hydrogen bonds formed between ATP and residues K52, Q103, D104, and M106. We expect that this study can provide a significant theoretical hint for designs of anticancer drugs targeting ERK2.