Force field independent metal parameters using a nonbonded dummy model

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.

Frontiers and Challenges of Computing ncRNAs Biogenesis, Function and Modulation

Non-coding RNAs (ncRNAs), generated from nonprotein coding DNA sequences, constitute 98-99% of the human genome. Non-coding RNAs encompass diverse functional classes, including microRNAs, small interfering RNAs, PIWI-interacting RNAs, small nuclear RNAs, small nucleolar RNAs, and long non-coding RNAs. With critical involvement in gene expression and regulation across various biological and physiopathological contexts, such as neuronal disorders, immune responses, cardiovascular diseases, and cancer, non-coding RNAs are emerging as disease biomarkers and therapeutic targets. In this review, after providing an overview of non-coding RNAs' role in cell homeostasis, we illustrate the potential and the challenges of state-of-the-art computational methods exploited to study non-coding RNAs biogenesis, function, and modulation. This can be done by directly targeting them with small molecules or by altering their expression by targeting the cellular engines underlying their biosynthesis. Drawing from applications, also taken from our work, we showcase the significance and role of computer simulations in uncovering fundamental facets of ncRNA mechanisms and modulation. This information may set the basis to advance gene modulation tools and therapeutic strategies to address unmet medical needs.

Computational investigation of copper-mediated conformational changes in α-synuclein dimer

We report molecular dynamics simulation of dimers of α-synuclein, the peptide closely associated with onset of Parkinson's disease, both as metal-free dimer and with inter-chain bridging provided by Cu(II) ions. Our investigation reveals that the presence of copper-induced inter-chain bridging not only stabilizes α-synuclein dimers, but also leads to enhanced β-sheet formation at critical regions within the N-terminal and NAC regions of the protein. These contacts are larger and longer-lived in the presence of copper, and as a result each peptide chain is more extended and less flexible than in the metal-free dimer. The persistence of these inter-peptide contacts underscores their significance in stabilising the dimers, potentially influencing the aggregation pathway. Moreover, the increased flexibility in the two termini, as well as the absence of persistent contacts in the metal-free dimer, correlates with the presence of amorphous aggregates. This phenomenon is known to mitigate fibrillation, while their absence in the metal-bound dimer suggests an increased propensity to form fibrils in the presence of copper ions.

DP/MM: A Hybrid Model for Zinc-Protein Interactions in Molecular Dynamics

Zinc-containing proteins are vital for many biological processes, yet accurately modeling them using classical force fields is hindered by complicated polarization and charge transfer effects. This study introduces DP/MM, a hybrid force field scheme that utilizes a deep potential model to correct the atomic forces of zinc ions and their coordinated atoms, elevating them from MM to QM levels of accuracy. Trained on the difference between MM and QM atomic forces across diverse zinc coordination groups, the DP/MM model faithfully reproduces structural characteristics of zinc coordination during simulations, such as the tetrahedral coordination of Cys4 and Cys3His1 groups. Furthermore, DP/MM allows water exchange in the zinc coordination environment. With its unique blend of accuracy, efficiency, flexibility, and transferability, DP/MM serves as a valuable tool for studying structures and dynamics of zinc-containing proteins and also represents a pioneering approach in the evolving landscape of machine learning potentials for molecular modeling.

Atomistic molecular simulations of Aβ-Zn conformational ensembles

The amyloid-forming Aβ peptide is able to interact with metal cations to form very stable complexes that influence fibril formation and contribute to the onset of Alzheimer's disease. Multiple structures of peptides derived from Aβ in complex with different metals have been resolved experimentally to provide an atomic-level description of the metal-protein interactions. However, Aβ is intrinsically disordered, and hence more amenable to an ensemble description. Molecular dynamics simulations can now reach the timescales needed to generate ensembles for these type of complexes. However, this requires accurate force fields both for the protein and the protein-metal interactions. Here we use state-of-the-art methods to generate force field parameters for the Zn(II) cations in a set of complexes formed by different Aβ variants and combine them with the Amber99SB*-ILDN optimized force field. Upon comparison of NMR experiments with the simulation results, further optimized with a Bayesian/Maximum entropy approach, we provide an accurate description of the molecular ensembles for most Aβ-metal complexes. We find that the resulting conformational ensembles are more heterogeneous than the NMR models deposited in the Protein Data Bank.

Attempting Well-Tempered Funnel Metadynamics Simulations for the Evaluation of the Binding Kinetics of Methionine Aminopeptidase-II Inhibitors

Understanding the unbinding kinetics of protein-ligand complexes is considered a significant approach for the design of ligands with desired specificity and safety. In recent years, enhanced sampling methods have emerged as effective tools for studying the unbinding kinetics of protein-ligand complexes at the atomistic level. MetAP-II is a target for the treatment of cancer for which not a single effective drug is available yet. The identification of the dissociation rate of ligands from the complexes often serves as a better predictor for in vivo efficacy than the ligands' binding affinity. Here, funnel-based restraint well-tempered metadynamics simulations were applied to predict the residence time of two ligands bound to MetAP-II, along with the ligand association and dissociation mechanism involving the identification of the binding hotspot during ligand egress. The ligand-egressing route revealed by metadynamics simulations also correlated with the identified pathways from the CAVER analysis and by the enhanced sampling simulation using PLUMED. Ligand 1 formed a strong H-bond interaction with GLU364 estimating a higher residence time of 28.22 ± 5.29 ns in contrast to ligand 2 with a residence time of 19.05 ± 3.58 ns, which easily dissociated from the binding pocket of MetAP-II. The results obtained from the simulations were consistent to reveal ligand 1 being superior to ligand 2; however, the experimental data related to residence time were close for both ligands, and no kinetic data were available for ligand 2. The current study could be considered the first attempt to apply an enhanced sampling method for the evaluation of the binding kinetics and thermodynamics of two different classes of ligands to a binuclear metalloprotein.

A sensor complements the steric gate when DNA polymerase ϵ discriminates ribonucleotides

The cellular imbalance between high concentrations of ribonucleotides (NTPs) and low concentrations of deoxyribonucleotides (dNTPs), is challenging for DNA polymerases when building DNA from dNTPs. It is currently believed that DNA polymerases discriminate against NTPs through a steric gate model involving a clash between a tyrosine and the 2'-hydroxyl of the ribonucleotide in the polymerase active site in B-family DNA polymerases. With the help of crystal structures of a B-family polymerase with a UTP or CTP in the active site, molecular dynamics simulations, biochemical assays and yeast genetics, we have identified a mechanism by which the finger domain of the polymerase sense NTPs in the polymerase active site. In contrast to the previously proposed polar filter, our experiments suggest that the amino acid residue in the finger domain senses ribonucleotides by steric hindrance. Furthermore, our results demonstrate that the steric gate in the palm domain and the sensor in the finger domain are both important when discriminating NTPs. Structural comparisons reveal that the sensor residue is conserved among B-family polymerases and we hypothesize that a sensor in the finger domain should be considered in all types of DNA polymerases.

Comparison of force fields to study the zinc-finger containing protein NPL4, a target for disulfiram in cancer therapy

Molecular dynamics (MD) simulations are a powerful approach to studying the structure and dynamics of proteins related to health and disease. Advances in the MD field allow modeling proteins with high accuracy. However, modeling metal ions and their interactions with proteins is still challenging. NPL4 is a zinc-binding protein and works as a cofactor for p97 to regulate protein homeostasis. NPL4 is of biomedical importance and has been proposed as the target of disulfiram, a drug recently repurposed for cancer treatment. Experimental studies proposed that the disulfiram metabolites, bis-(diethyldithiocarbamate)‑copper and cupric ions, induce NPL4 misfolding and aggregation. However, the molecular details of their interactions with NPL4 and consequent structural effects are still elusive. Here, biomolecular simulations can help to shed light on the related structural details. To apply MD simulations to NPL4 and its interaction with copper the first important step is identifying a suitable force field to describe the protein in its zinc-bound states. We examined different sets of non-bonded parameters because we want to study the misfolding mechanism and cannot rule out that the zinc may detach from the protein during the process and copper replaces it. We investigated the force-field ability to model the coordination geometry of the metal ions by comparing the results from MD simulations with optimized geometries from quantum mechanics (QM) calculations using model systems of NPL4. Furthermore, we investigated the performance of a force field including bonded parameters to treat copper ions in NPL4 that we obtained based on QM calculations.

MDiGest: A Python package for describing allostery from molecular dynamics simulations

Many biological processes are regulated by allosteric mechanisms that communicate with distant sites in the protein responsible for functionality. The binding of a small molecule at an allosteric site typically induces conformational changes that propagate through the protein along allosteric pathways regulating enzymatic activity. Elucidating those communication pathways from allosteric sites to orthosteric sites is, therefore, essential to gain insights into biochemical processes. Targeting the allosteric pathways by mutagenesis can allow the engineering of proteins with desired functions. Furthermore, binding small molecule modulators along the allosteric pathways is a viable approach to target reactions using allosteric inhibitors/activators with temporal and spatial selectivity. Methods based on network theory can elucidate protein communication networks through the analysis of pairwise correlations observed in molecular dynamics (MD) simulations using molecular descriptors that serve as proxies for allosteric information. Typically, single atomic descriptors such as α-carbon displacements are used as proxies for allosteric information. Therefore, allosteric networks are based on correlations revealed by that descriptor. Here, we introduce a Python software package that provides a comprehensive toolkit for studying allostery from MD simulations of biochemical systems. MDiGest offers the ability to describe protein dynamics by combining different approaches, such as correlations of atomic displacements or dihedral angles, as well as a novel approach based on the correlation of Kabsch-Sander electrostatic couplings. MDiGest allows for comparisons of networks and community structures that capture physical information relevant to allostery. Multiple complementary tools for studying essential dynamics include principal component analysis, root mean square fluctuation, as well as secondary structure-based analyses.

Lenalidomide Stabilizes Protein-Protein Complexes by Turning Labile Intermolecular H-Bonds into Robust Interactions

Targeted protein degradation is a promising therapeutic strategy, spearheaded by the anti-myeloma drugs lenalidomide and pomalidomide. These drugs stabilize very efficiently the complex between the E3 ligase Cereblon (CRBN) and several non-native client proteins (neo-substrates), including the transcription factors Ikaros and Aiolos and the enzyme Caseine Kinase 1α (CK1α,), resulting in their degradation. Although the structures for these complexes have been determined, there are no evident interactions that can account for the high efficiency of formation of the ternary complex. We show that lenalidomide's stabilization of the CRBN-CK1α complex is largely due to hydrophobic shielding of intermolecular hydrogen bonds. We also find a quantitative relationship between hydrogen bond robustness and binding affinities of the ternary complexes. These results pave the way to further understand cooperativity effects in drug-induced protein-protein complexes and could help in the design of improved molecular glues and more efficient protein degraders.

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.

Modeling of minimal systems based on ATP-Zn coordination for chemically fueled self-assembly

Following nature's example, there is currently strong interest in using adenosine 5'-triphosphate (ATP) as a fuel for the self-assembly of functional materials with transient/non-equilibrium behaviours. These hold great promise for applications, e.g. in catalysis and drug delivery. In a recent seminal work [Maiti et al., Nat. Chem., 2016, 8, 725], binding of ATP to the metallosurfactant zinc hexadecyl-1,4,7-triazacyclononane ([ZnC₁₆ TACN]²⁺) was exploited to produce ATP-fueled transient vesicles. Crucial to the complex formation is the ability of ATP to bind to the metal ion. As a first step to unveil the key elements underlying this process, we investigate the interaction of ATP with Zn²⁺ and with methyl-1,4,7-triazacyclononane ([ZnCH₃ TACN]²⁺), using all-atom molecular dynamics simulations. The free energy landscape of the complex formation is sampled using well-tempered metadynamics with three collective variables, corresponding to the coordination numbers of Zn²⁺ with the oxygen atoms of the three phosphate groups. We find that the structure of the ternary complex is controlled by direct triphosphate coordination to zinc, with a minor role played by the interactions between ATP and CH₃ TACN which, however, may be important for the build-up of supramolecular assemblies.

Allosteric rescue of catalytically impaired ATP phosphoribosyltransferase variants links protein dynamics to active-site electrostatic preorganisation

ATP phosphoribosyltransferase catalyses the first step of histidine biosynthesis and is controlled via a complex allosteric mechanism where the regulatory protein HisZ enhances catalysis by the catalytic protein HisG_S while mediating allosteric inhibition by histidine. Activation by HisZ was proposed to position HisG_S Arg56 to stabilise departure of the pyrophosphate leaving group. Here we report active-site mutants of HisG_S with impaired reaction chemistry which can be allosterically restored by HisZ despite the HisZ:HisG_S interface lying ~20 Å away from the active site. MD simulations indicate HisZ binding constrains the dynamics of HisG_S to favour a preorganised active site where both Arg56 and Arg32 are poised to stabilise leaving-group departure in WT-HisG_S. In the Arg56Ala-HisG_S mutant, HisZ modulates Arg32 dynamics so that it can partially compensate for the absence of Arg56. These results illustrate how remote protein-protein interactions translate into catalytic resilience by restoring damaged electrostatic preorganisation at the active site.

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.

A Polarizable Cationic Dummy Metal Ion Model

A novel locally polarizable multisite model based on the original cation dummy atom (CDA) model is described for molecular dynamics simulations of ions in condensed phases. Polarization effects are introduced by the electronegativity equalization model (EEM) method where charges on the metal ion and its dummy atoms can fluctuate to respond to the environment. This model includes explicit polarization and ion-induced interactions and can be coupled with nonpolarizable or polarizable water models, making it more transferable to simpler force fields. This approach allows us to enhance the original fixed charge CDA model where the charge distribution cannot adapt to the local solvent structure. To illustrate the new CDApol model, we examined properties of the Zn2+, Al3+, and Zr4+ ions in aqueous solution. The polarizable model and Lennard-Jones parameters were refined for octahedrally coordinated Zn2+, Al3+, and Zr4+ CDAs to reproduce thermodynamic and geometrical properties. Using this locally polarizable model, we were able to obtain the experimental hydration free energy, ion-oxygen distance, and coordination number coupled with the standard 12-6 Lennard-Jones model. This model can be used in myriad additional applications where local polarization and charge transfer is important.

Metal Ion Binding in Wild-Type and Mutated Frataxin: A Stability Study

This work studies the stability of wild-type frataxin and some of its variants found in cancer tissues upon Co2+ binding. Although the physiologically involved metal ion in the frataxin enzymatic activity is Fe2+, as it is customarily done, Co2+ is most often used in experiments because Fe2+ is extremely unstable owing to the fast oxidation reaction Fe2+ → Fe3+. Protein stability is monitored following the conformational changes induced by Co2+ binding as measured by circular dichroism, fluorescence spectroscopy, and melting temperature measurements. The stability ranking among the wild-type frataxin and its variants obtained in this way is confirmed by a detailed comparative analysis of the XAS spectra of the metal-protein complex at the Co K-edge. In particular, a fit to the EXAFS region of the spectrum allows positively identifying the frataxin acidic ridge as the most likely location of the metal-binding sites. Furthermore, we can explain the surprising feature emerging from a detailed analysis of the XANES region of the spectrum, showing that the longer 81-210 frataxin fragment has a smaller propensity for Co2+ binding than the shorter 90-210 one. This fact is explained by the peculiar role of the N-terminal disordered tail in modulating the protein ability to interact with the metal.

Structure-Based Discovery of Lipoteichoic Acid Synthase Inhibitors

Lipoteichoic acid synthase (LtaS) is a key enzyme for the cell wall biosynthesis of Gram-positive bacteria. Gram-positive bacteria that lack lipoteichoic acid (LTA) exhibit impaired cell division and growth defects. Thus, LtaS appears to be an attractive antimicrobial target. The pharmacology around LtaS remains largely unexplored with only two small-molecule LtaS inhibitors reported, namely "compound 1771" and the Congo red dye. Structure-based drug discovery efforts against LtaS remain unattempted due to the lack of an inhibitor-bound structure of LtaS. To address this, we combined the use of a molecular docking technique with molecular dynamics (MD) simulations to model a plausible binding mode of compound 1771 to the extracellular catalytic domain of LtaS (eLtaS). The model was validated using alanine mutagenesis studies combined with isothermal titration calorimetry. Additionally, lead optimization driven by our computational model resulted in an improved version of compound 1771, namely, compound 4 which showed greater affinity for binding to eLtaS than compound 1771 in biophysical assays. Compound 4 reduced LTA production in S. aureus dose-dependently, induced aberrant morphology as seen for LTA-deficient bacteria, and significantly reduced bacteria titers in the lung of mice infected with S. aureus. Analysis of our MD simulation trajectories revealed the possible formation of a transient cryptic pocket in eLtaS. Virtual screening (VS) against the cryptic pocket led to the identification of a new class of inhibitors that could potentiate β-lactams against methicillin-resistant S. aureus. Our overall workflow and data should encourage further drug design campaign against LtaS. Finally, our work reinforces the importance of considering protein conformational flexibility to a successful VS endeavor.

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.

A multiscale approach to predict the binding mode of metallo beta-lactamase inhibitors

Antibiotic resistance is a major threat to global public health. β-lactamases, which catalyze breakdown of β-lactam antibiotics, are a principal cause. Metallo β-lactamases (MBLs) represent a particular challenge because they hydrolyze almost all β-lactams and to date no MBL inhibitor has been approved for clinical use. Molecular simulations can aid drug discovery, for example, predicting inhibitor complexes, but empirical molecular mechanics (MM) methods often perform poorly for metalloproteins. Here we present a multiscale approach to model thiol inhibitor binding to IMP-1, a clinically important MBL containing two catalytic zinc ions, and predict the binding mode of a 2-mercaptomethyl thiazolidine (MMTZ) inhibitor. Inhibitors were first docked into the IMP-1 active site, testing different docking programs and scoring functions on multiple crystal structures. Complexes were then subjected to molecular dynamics (MD) simulations and subsequently refined through QM/MM optimization with a density functional theory (DFT) method, B3LYP/6-31G(d), increasing the accuracy of the method with successive steps. This workflow was tested on two IMP-1:MMTZ complexes, for which it reproduced crystallographically observed binding, and applied to predict the binding mode of a third MMTZ inhibitor for which a complex structure was crystallographically intractable. We also tested a 12-6-4 nonbonded interaction model in MD simulations and optimization with a SCC-DFTB QM/MM approach. The results show the limitations of empirical models for treating these systems and indicate the need for higher level calculations, for example, DFT/MM, for reliable structural predictions. This study demonstrates a reliable computational pipeline that can be applied to inhibitor design for MBLs and other zinc-metalloenzyme systems.

Hydration modulates oxygen channel residues for oxygenation of cysteine dioxygenase: Perspectives from molecular dynamics simulations

Cysteine dioxygenase (CDO) regulates the concentration of l-cysteine substrate by its oxidation in the body to prevent different diseases, including neurodegenerative and autoimmune diseases. CDO catalyzes the oxidation of thiol group of l-cysteine to l-cysteine sulfinic acid using molecular oxygen. In this study, molecular dynamics simulations were applied to ligand-free CDO, cysteine-bound CDO, and oxygen-bound CDO-cysteine complex which were primarily subjected to the evaluation of their structural and dynamical properties. The simulation data provided significant information not only on the conformational changes of the enzyme after its ligation but also on the co-ligation by sequential binding of l-cysteine and molecular oxygen. It was found that the ligation and co-ligation perturbed the active site region as well as the overall protein dynamics which were analyzed in terms of root mean square deviation, root mean square fluctuation and dynamic cross correlation matrices as well as principal component analysis. Furthermore, oxygen transport pathways were successfully explored by taking various tunnel clusters into account and one of those clusters was given preference based on the throughput value. The bottleneck formed by different amino acid residues was examined to figure out their role in the oxygenation process of the enzyme. The residues forming the tunnel's bottleneck and their dynamics mediated by water molecules were further investigated using radial distribution functions which gave insights into the hydration behavior of these residues. The findings based on the hydration behavior in turn served to explore the water-mediated dynamics of these residues in the modulation of the pathway, including tunnel gating for the oxygen entry and diffusion to the active site, which is essential for the CDO's catalytic function.

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.

Modeling Kinetics and Thermodynamics of Guest Encapsulation into the [M4L6]12- Supramolecular Organometallic Cage

The encapsulation of molecular guests into supramolecular hosts is a complex molecular recognition process in which the guest displaces the solvent from the host cavity, while the host deforms to let the guest in. An atomistic description of the association would provide valuable insights on the physicochemical properties that guide it. This understanding may be used to design novel host assemblies with improved properties (i.e., affinities) toward a given class of guests. Molecular simulations may be conveniently used to model the association processes. It is thus of interest to establish efficient protocols to trace the encapsulation process and to predict the associated magnitudes ΔG_bind and ΔG_bind^⧧. Here, we report the calculation of the Gibbs energy barrier and Gibbs binding energy by means of explicit solvent molecular simulations for the [Ga₄L₆]^12– metallocage encapsulating a series of cationic molecules. The ΔG_bind^⧧ for encapsulation was estimated by means of umbrella sampling simulations. The steps involved were identified, including ion-pair formation and naphthalene rotation (from L ligands of the metallocage) during the guest’s entrance. The ΔG_bind values were computed using the attach–pull–release method. The results reveal the sensitivity of the estimates on the force field parameters, in particular on atomic charges, showing that higher accuracy is obtained when charges are derived from implicit solvent quantum chemical calculations. Correlation analysis identified some indicators for the binding affinity trends. All computed magnitudes are in very good agreement with experimental observations. This work provides, on one side, a benchmarked way to computationally model a highly charged metallocage encapsulation process. This includes a nonstandard parameterization and charge derivation procedure. On the other hand, it gives specific mechanistic information on the binding processes of [Ga₄L₆]^12– at the molecular level where key motions are depicted. Taken together, the study provides an interesting option for the future design of metal–organic cages.

Probing the Suitability of Different Ca2+ Parameters for Long Simulations of Diisopropyl Fluorophosphatase

Organophosphate hydrolases are promising as potential biotherapeutic agents to treat poisoning with pesticides or nerve gases. However, these enzymes often need to be further engineered in order to become useful in practice. One example of such enhancement is the alteration of enantioselectivity of diisopropyl fluorophosphatase (DFPase). Molecular modeling techniques offer a unique opportunity to address this task rationally by providing a physical description of the substrate-binding process. However, DFPase is a metalloenzyme, and correct modeling of metal cations is a challenging task generally coming with a tradeoff between simulation speed and accuracy. Here, we probe several molecular mechanical parameter combinations for their ability to empower long simulations needed to achieve a quantitative description of substrate binding. We demonstrate that a combination of the Amber19sb force field with the recently developed 12-6 Ca2+ models allows us to both correctly model DFPase and obtain new insights into the DFP binding process.

Coarse-Grained Modeling and Molecular Dynamics Simulations of Ca2+-Calmodulin

Calmodulin (CaM) is a calcium-binding protein that transduces signals to downstream proteins through target binding upon calcium binding in a time-dependent manner. Understanding the target binding process that tunes CaM’s affinity for the calcium ions (Ca²⁺), or vice versa, may provide insight into how Ca²⁺-CaM selects its target binding proteins. However, modeling of Ca²⁺-CaM in molecular simulations is challenging because of the gross structural changes in its central linker regions while the two lobes are relatively rigid due to tight binding of the Ca²⁺ to the calcium-binding loops where the loop forms a pentagonal bipyramidal coordination geometry with Ca²⁺. This feature that underlies the reciprocal relation between Ca²⁺ binding and target binding of CaM, however, has yet to be considered in the structural modeling. Here, we presented a coarse-grained model based on the Associative memory, Water mediated, Structure, and Energy Model (AWSEM) protein force field, to investigate the salient features of CaM. Particularly, we optimized the force field of CaM and that of Ca²⁺ ions by using its coordination chemistry in the calcium-binding loops to match with experimental observations. We presented a “community model” of CaM that is capable of sampling various conformations of CaM, incorporating various calcium-binding states, and carrying the memory of binding with various targets, which sets the foundation of the reciprocal relation of target binding and Ca²⁺ binding in future studies.

A novel mechanism of enhanced transcription activity and fidelity for influenza A viral RNA-dependent RNA polymerase

During RNA elongation, the influenza A viral (IAV) RNA-dependent RNA polymerase (RdRp) residues in the active site interact with the triphosphate moiety of nucleoside triphosphate (NTP) for catalysis. The molecular mechanisms by which they control the rate and fidelity of NTP incorporation remain elusive. Here, we demonstrated through enzymology, virology and computational approaches that the R239 and K235 in the PB1 subunit of RdRp are critical to controlling the activity and fidelity of transcription. Contrary to common beliefs that high-fidelity RdRp variants exert a slower incorporation rate, we discovered a first-of-its-kind, single lysine-to-arginine mutation on K235 exhibited enhanced fidelity and activity compared with wild-type. In particular, we employed a single-turnover NTP incorporation assay for the first time on IAV RdRp to show that K235R mutant RdRp possessed a 1.9-fold increase in the transcription activity of the cognate NTP and a 4.6-fold increase in fidelity compared to wild-type. Our all-atom molecular dynamics simulations further elucidated that the higher activity is attributed to the shorter distance between K235R and the triphosphate moiety of NTP compared with wild-type. These results provide novel insights into NTP incorporation and fidelity control mechanisms, which lay the foundation for the rational design of IAV vaccine and antiviral targets.

Unravelling the inhibitory zinc ion binding site and the metal exchange mechanism in human DPP III

Dipeptidyl peptidase III (DPP III), a zinc-dependent exopeptidase, is widely distributed in organisms and present in almost all human tissues. In addition to its involvement in protein catabolism, it plays a role in oxidative stress and blood pressure regulation, and there is evidence of its involvement in pain modulation. Excess zinc ions have been found to inhibit its hydrolytic activity, but the binding affinity, binding site geometry, and mechanism of inhibitory activity have been unknown. Using several different computational approaches, we determined the inhibitory zinc ion binding site, its coordination and relative binding affinity. During some simulations the translocation of the zinc ion from the inhibitory to the catalytic binding site was observed, accompanied by movement of the catalytic zinc ion toward the exit of the substrate binding site. The traced behavior suggests an associative type of metal ion exchange, in which the formation of the ternary complex between enzyme and two metal ions precedes the exit of the catalytic metal ion. Differently from our previous findings that binding of a peptide induces partial opening of hDPP III, the globularity of the protein did not change in MD simulations of the hermorphin-like peptide bound to hDPP III with two zinc ions. However, the entrance to the interdomain cleft widens during Zn diffusion into the protein and was found to be the highest energy barrier in the process of metal translocation from the solvent to the active site. Finally, we discuss why excess zinc reduces enzyme activity.

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.

Zinc binding alters the conformational dynamics and drives the transport cycle of the cation diffusion facilitator YiiP

YiiP is a secondary transporter that couples Zn2+ transport to the proton motive force. Structural studies of YiiP from prokaryotes and Znt8 from humans have revealed three different Zn2+ sites and a conserved homodimeric architecture. These structures define the inward-facing and outward-facing states that characterize the archetypal alternating access mechanism of transport. To study the effects of Zn2+ binding on the conformational transition, we use cryo-EM together with molecular dynamics simulation to compare structures of YiiP from Shewanella oneidensis in the presence and absence of Zn2+. To enable single-particle cryo-EM, we used a phage-display library to develop a Fab antibody fragment with high affinity for YiiP, thus producing a YiiP/Fab complex. To perform MD simulations, we developed a nonbonded dummy model for Zn2+ and validated its performance with known Zn2+-binding proteins. Using these tools, we find that, in the presence of Zn2+, YiiP adopts an inward-facing conformation consistent with that previously seen in tubular crystals. After removal of Zn2+ with high-affinity chelators, YiiP exhibits enhanced flexibility and adopts a novel conformation that appears to be intermediate between inward-facing and outward-facing states. This conformation involves closure of a hydrophobic gate that has been postulated to control access to the primary transport site. Comparison of several independent cryo-EM maps suggests that the transition from the inward-facing state is controlled by occupancy of a secondary Zn2+ site at the cytoplasmic membrane interface. This work enhances our understanding of individual Zn2+ binding sites and their role in the conformational dynamics that govern the transport cycle.

The mechanism of action of T-705 as a unique delayed chain terminator on influenza viral polymerase transcription

Favipiravir (T-705) has been developed as a potent anti-influenza drug and exhibited a strong inhibition effect against a broad spectrum of RNA viruses. Its active form, ribofuranosyl-triphosphate (T-705-RTP), functions as a competitive substrate for the RNA-dependent RNA polymerase (RdRp) of the influenza A virus (IAV). However, the exact inhibitory mechanisms of T-705 remain elusive and subject to a long-standing debate. Although T-705 has been proposed to inhibit transcription by acting as a chain terminator, it is also paradoxically suggested to be a mutagen towards IAV RdRp by inducing mutations due to its ambiguous base pairing of C and U. Here, we combined biochemical assay with molecular dynamics (MD) simulations to elucidate the molecular mechanism underlying the inhibitory functions exerted by T-705 in IAV RdRp. Our in vitro transcription assay illustrated that IAV RdRp could recognize T-705 as a purine analogue and incorporate it into the nascent RNA strand. Incorporating a single T-705 is incapable of inhibiting transcription as extra natural nucleotides can be progressively added. However, when two consecutive T-705 are incorporated, viral transcription is completely terminated. MD simulations reveal that the sequential appearance of two T-705 in the nascent strand destabilizes the active site and disrupts the base stacking of the nascent RNA. Altogether, our results provide a plausible explanation for the inhibitory roles of T-705 targeting IAV RdRp by integrating the computational and experimental methods. Our study also offers a comprehensive platform to investigate the inhibition effect of antivirals and a novel explanation for the designing of anti-flu drugs.

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.

Gonçalves Vianez JL, Santana da Costa K, de Molfetta FA, Nahum Alves C. Analysis of Kojic Acid Derivatives as Competitive Inhibitors of Tyrosinase: A Molecular Modeling Approach

Tyrosinases belong to the functional copper-containing proteins family, and their structure contains two copper atoms, in the active site, which are coordinated by three histidine residues. The biosynthesis of melanin in melanocytes has two stages depending on the actions of the natural substrates L-DOPA and L-tyrosine. The dysregulation of tyrosinase is involved in skin cancer initiation. In the present study, using molecular modeling tools, we analyzed the inhibition activity of tyrosinase activity using kojic acid (KA) derivatives designed from aromatic aldehydes and malononitrile. All derivatives showed conformational affinity to the enzyme active site, and a favorable distance to chelate the copper ion, which is essential for enzyme function. Molecular dynamics simulations revealed that the derivatives formed promising complexes, presenting stable conformations with deviations between 0.2 and 0.35 Å. In addition, the investigated KA derivatives showed favorable binding free energies. The most stable KA derivatives showed the following binding free energies: -17.65 kcal mol-1 (D6), -18.07 kcal mol-1 (D2), -18.13 (D5) kcal mol-1, and -10.31 kcal mol-1 (D4). Our results suggest that these derivatives could be potent competitive inhibitors of the natural substrates of L-DOPA (-12.84 kcal mol-1) and L-tyrosine (-9.04 kcal mol-1) in melanogenesis.

Computing Metal-Binding Proteins for Therapeutic Benefit

Over one third of biomolecules rely on metal ions to exert their cellular functions. Metal ions can play a structural role by stabilizing the structure of biomolecules, a functional role by promoting a wide variety of biochemical reactions, and a regulatory role by acting as messengers upon binding to proteins regulating cellular metal-homeostasis. These diverse roles in biology ascribe critical implications to metal-binding proteins in the onset of many diseases. Hence, it is of utmost importance to exhaustively unlock the different mechanistic facets of metal-binding proteins and to harness this knowledge to rationally devise novel therapeutic strategies to prevent or cure pathological states associated with metal-dependent cellular dysfunctions. In this compendium, we illustrate how the use of a computational arsenal based on docking, classical, and quantum-classical molecular dynamics simulations can contribute to extricate the minutiae of the catalytic, transport, and inhibition mechanisms of metal-binding proteins at the atomic level. This knowledge represents a fertile ground and an essential prerequisite for selectively targeting metal-binding proteins with small-molecule inhibitors aiming to (i) abrogate deregulated metal-dependent (mis)functions or (ii) leverage metal-dyshomeostasis to selectively trigger harmful cells death.

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

Commonly seen in rare-earth chemistry and materials science, highly charged metal ions play key roles in many chemical processes. Computer simulations have become an important tool for scientific research nowadays. Meaningful simulations require reliable parameters. In the present work, we parametrized 18 M(III) and 6 M(IV) metal ions for four new water models (OPC3, OPC, TIP3P-FB, TIP4P-FB) in conjunction with each of the 12-6 and 12-6-4 nonbonded models. Similar to what was observed previously, issues with the 12-6 model can be fixed by using the 12-6-4 model. Moreover, the four new water models showed comparable performance or considerable improvement over the previous water models (TIP3P, SPC/E, and TIP4P_EW) in the same category (3-point or 4-point water models, respectively). Finally, we reported a study of a metalloprotein system demonstrating the capability of the 12-6-4 model to model metalloproteins. The reported parameters will facilitate accurate simulations of highly charged metal ions in aqueous solution.

A practical guide to biologically relevant molecular simulations with charge scaling for electronic polarization

Molecular simulations can elucidate atomistic-level mechanisms of key biological processes, which are often hardly accessible to experiment. However, the results of the simulations can only be as trustworthy as the underlying simulation model. In many of these processes, interactions between charged moieties play a critical role. Current empirical force fields tend to overestimate such interactions, often in a dramatic way, when polyvalent ions are involved. The source of this shortcoming is the missing electronic polarization in these models. Given the importance of such biomolecular systems, there is great interest in fixing this deficiency in a computationally inexpensive way without employing explicitly polarizable force fields. Here, we review the electronic continuum correction approach, which accounts for electronic polarization in a mean-field way, focusing on its charge scaling variant. We show that by pragmatically scaling only the charged molecular groups, we qualitatively improve the charge-charge interactions without extra computational costs and benefit from decades of force field development on biomolecular force fields.

Parameterization of Monovalent Ions for the OPC3, OPC, TIP3P-FB, and TIP4P-FB Water Models

Monovalent ions play significant roles in various biological and material systems. Recently, four new water models (OPC3, OPC, TIP3P-FB, and TIP4P-FB), with significantly improved descriptions of condensed phase water, have been developed. The pairwise interaction between the metal ion and water necessitates the development of ion parameters specifically for these water models. Herein, we parameterized the 12-6 and the 12-6-4 nonbonded models for 12 monovalent ions with the respective four new water models. These monovalent ions contain eight cations including alkali metal ions (Li⁺, Na⁺, K⁺, Rb⁺, Cs⁺), transition-metal ions (Cu⁺ and Ag⁺), and Tl⁺ from the boron family, along with four halide anions (F^-, Cl^-, Br^-, I^-). Our parameters were designed to reproduce the target hydration free energies (the 12-6 hydration free energy (HFE) set), the ion-oxygen distances (the 12-6 ion-oxygen distance (IOD) set), or both of them (the 12-6-4 set). The 12-6-4 parameter set provides highly accurate structural features overcoming the limitations of the routinely used 12-6 nonbonded model for ions. Specifically, we note that the 12-6-4 parameter set is able to reproduce experimental hydration free energies within 1 kcal/mol and experimental ion-oxygen distances within 0.01 Å simultaneously. We further reproduced the experimentally determined activity derivatives for salt solutions, validating the ion parameters for simulations of ion pairs. The improved performance of the present water models over our previous parameter sets for the TIP3P, TIP4P, and SPC/E water models (Li, P. et al J. Chem. Theory Comput. 2015 11 1645 1657) highlights the importance of the choice of water model in conjunction with the metal ion parameter set.

Evaluating the Performance of a Non-Bonded Cu2+ Model Including Jahn-Teller Effect into the Binding of Tyrosinase Inhibitors

Tyrosinase (TYR) is a metalloenzyme classified as a type-3 copper protein, which is involved in the synthesis of melanin through a catalytic process beginning with the conversion of the amino acid l-Tyrosine (l-Tyr) to l-3,4-dihydroxyphenylalanine (l-DOPA). It plays an important role in the mechanism of melanogenesis in various organisms including mammals, plants, and fungi. Herein, we used a combination of computational molecular modeling techniques including molecular dynamic (MD) simulations and the linear interaction energy (LIE) model to evaluate the binding free energy of a set of analogs of kojic acid (KA) in complex with TYR. For the MD simulations, we used a dummy model including the description of the Jahn–Teller effect for Cu²⁺ ions in the active site of this enzyme. Our results show that the LIE model predicts the TYR binding affinities of the inhibitor in close agreement to experimental results. Overall, we demonstrate that the classical model provides a suitable description of the main interactions between analogs of KA and Cu²⁺ ions in the active site of TYR.

Upgraded AMBER Force Field for Zinc-Binding Residues and Ligands for Predicting Structural Properties and Binding Affinities in Zinc-Proteins

We developed a novel force field in the context of AMBER parameterization for glutamate and aspartate zinc(II)-binding residues. The interaction between the zinc ion and the coordinating atoms is represented by a spherical nonbonded parameterization. The polarization effect due to the zinc ion has been taken into account by redefining the atomic charges on the residues through accurate quantum mechanical calculations. The new zinc-binding ASP and GLU residues, along with the CYS and HIS zinc-binding residues, parameterized in a recent work [Macchiagodena M.;J. Chem. Inf. Model.2019, 59, 3803-3816], allow users to reliably simulate 96% of the Zn-proteins available in the Protein Data Bank. The upgraded force field for zinc(II)-bound residues has been tested performing molecular dynamics simulations with an explicit solvent and comparing the structural information with experimental data for five different proteins binding zinc(II) with GLU, ASP, HIS, and CYS. We further validated our approach by evaluating the binding free energy of (R)-2-benzyl-3-nitropropanoic acid to carboxypeptidase A using a recently developed nonequilibrium alchemical method. We demonstrated that in this setting it is crucial to take into account polarization effects also on the metal-bound inhibitor.

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.

A kinase bioscavenger provides antibiotic resistance by extremely tight substrate binding

Microbial communities are self-controlled by repertoires of lethal agents, the antibiotics. In their turn, these antibiotics are regulated by bioscavengers that are selected in the course of evolution. Kinase-mediated phosphorylation represents one of the general strategies for the emergence of antibiotic resistance. A new subfamily of AmiN-like kinases, isolated from the Siberian bear microbiome, inactivates antibiotic amicoumacin by phosphorylation. The nanomolar substrate affinity defines AmiN as a phosphotransferase with a unique catalytic efficiency proximal to the diffusion limit. Crystallographic analysis and multiscale simulations revealed a catalytically perfect mechanism providing phosphorylation exclusively in the case of a closed active site that counteracts substrate promiscuity. AmiN kinase is a member of the previously unknown subfamily representing the first evidence of a specialized phosphotransferase bioscavenger.

Molecular Modeling for Artificial Metalloenzyme Design and Optimization

Artificial metalloenzymes (ArMs) are obtained by inserting homogeneous catalysts into biological scaffolds and are among the most promising strategies in the quest for new-to-nature biocatalysts. The quality of their design strongly depends on how three partners interact: the biological host, the "artificial cofactor," and the substrate. However, structural characterization of functional artificial metalloenzymes by X-ray or NMR is often partial, elusive, or absent. How the cofactor binds to the protein, how the receptor reorganizes upon the binding of the cofactor and the substrate, and which are the binding mode(s) of the substrate for the reaction to proceed are key questions that are frequently unresolved yet crucial for ArM design. Such questions may eventually be solved by molecular modeling but require a step change beyond the current state-of-the-art methodologies.Here, we summarize our efforts in the study of ArMs, presenting both the development of computational strategies and their application. We first focus on our integrative computational framework that incorporates a variety of methods such as protein-ligand docking, classical molecular dynamics (MD), and pure quantum mechanical (QM) methods, which, when properly combined, are able to depict questions that range from host-cofactor binding predictions to simulations of entire catalytic mechanisms. We also pay particular attention to the protein-ligand docking strategies that we have developed to accurately predict the binding of transition metal-containing molecules to proteins. While this aspect is fundamental to many bioinorganic fields beyond ArMs, it has been disregarded from the molecular modeling landscape until very recently.Next we describe how to apply this computational framework to particular ArMs including systems previously characterized experimentally as well as others where computation served to guide the design. We start with the prediction of the interactions between homogeneous catalysts and biological hosts. Protein-ligand docking is pivotal at that stage, but it needs to be combined with QM/MM or MD approaches when the binding of the cofactor implies significant conformational changes of the protein or involve changes of the electronic state of the metal.Then, we summarize molecular modeling studies aimed at identifying cofactor-substrate arrangements inside the ArM active pocket that are consistent with its reactivity. These calculations stand on "Theozyme"-like dockings, MD-refined or not, which provide molecular rationale of the catalytic profiles of the artificial systems.In the third section, we present case studies to decode the entire catalytic mechanism of two ArMs: (1) an iridium based asymmetric transfer hydrogenase obtained by insertion of Noyori's catalyst into streptavidin and (2) a metallohydrolase achieved by including a receptor. Transition states, second coordination sphere effects, as well as motions of the cofactors are identified as drivers of the enantiomeric profiles.Finally, we report computer-aided designs of ArMs to guide experiments toward chemical and mutational changes that improve their activity and/or enantioselective profiles and expand toward future directions.

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.

The Ca2+ permeation mechanism of the ryanodine receptor revealed by a multi-site ion model

Ryanodine receptors (RyR) are ion channels responsible for the release of Ca²⁺ from the sarco/endoplasmic reticulum and play a crucial role in the precise control of Ca²⁺ concentration in the cytosol. The detailed permeation mechanism of Ca²⁺ through RyR is still elusive. By using molecular dynamics simulations with a specially designed Ca²⁺ model, we show that multiple Ca²⁺ ions accumulate in the upper selectivity filter of RyR1, but only one Ca²⁺ can occupy and translocate in the narrow pore at a time, assisted by electrostatic repulsion from the Ca²⁺ within the upper selectivity filter. The Ca²⁺ is nearly fully hydrated with the first solvation shell intact during the whole permeation process. These results suggest a remote knock-on permeation mechanism and one-at-a-time occupation pattern for the hydrated Ca²⁺ within the narrow pore, uncovering the basis underlying the high permeability and low selectivity of the RyR channels.

AMBER Force Field Parameters for Cobalt-Containing Biological Systems: A Systematic Derivation Study

In the present work, the parameterization of a set of cobalt-containing systems has been performed to create a comprehensive library for bonding parameters of biological Co-containing systems. A standard process for the extraction and validation of parameters was employed, which could be used to create force field parameters for the other metal-containing systems. All protein data banks were searched to extract common chemical groups in bonding with cobalt, and finally, 16 structures were designed to represent the binding model of the chemical moieties with cobalt. The Hessian matrix of each structure was computed at the B3LYP/6-311++G(2d,2p) level of theory and the Seminario method was employed to compute cobalt bond stretching and angle bending parameters. Validation of the derived parameters was performed using structural minimization and molecular dynamics (MD) simulations of four models. Further validation was performed using an extensive MD simulation on carbonic anhydrase II as a common cobalt-containing metalloprotein. The results demonstrated that among models, the bonded model in combination with the RESP charges can produce the most reliable and accurate structural conformations for the metal site of cobalt-containing systems.