Force fields including charge transfer and local polarization effects: Application to proteins containing multi/heavy metal ions

System-Specific Parameter Optimization for Nonpolarizable and Polarizable Force Fields

The accuracy of classical force fields (FFs) has been shown to be limited for the simulation of cation-protein systems despite their importance in understanding the processes of life. Improvements can result from optimizing the parameters of classical FFs or by extending the FF formulation by terms describing charge transfer (CT) and polarization (POL) effects. In this work, we introduce our implementation of the CTPOL model in OpenMM, which extends the classical additive FF formula by adding CT and POL. Furthermore, we present an open-source parametrization tool, called FFAFFURR, that enables the (system-specific) parametrization of OPLS-AA and CTPOL models. The performance of our workflow was evaluated by its ability to reproduce quantum chemistry energies and by molecular dynamics simulations of a zinc-finger protein.

Protein Ca2+-Sites Prone to Sr2+ Substitution: Implications for Strontium Therapy

Strontium (Sr), an alkali metal with properties similar to calcium, in the form of soluble salts is used to treat osteoporosis. Despite the information accumulated on the role of Sr²⁺ as a Ca²⁺ mimetic in biology and medicine, there is no systematic study of how the outcome of the competition between the two dications depends on the physicochemical properties of (i) the metal ions, (ii) the first- and second-shell ligands, and (iii) the protein matrix. Specifically, the key features of a Ca²⁺-binding protein that enable Sr²⁺ to displace Ca²⁺ remain unclear. To address this, we studied the competition between Ca²⁺ and Sr²⁺ in protein Ca²⁺-binding sites using density functional theory combined with the polarizable continuum model. Our findings indicate that Ca²⁺-sites with multiple strong charge-donating protein ligands, including one or more bidentately bound Asp^-/Glu^- that are relatively buried and rigid are protected against Sr²⁺ attack. On the other hand, Ca²⁺-sites crowded with multiple protein ligands may be prone to Sr²⁺ displacement if they are solvent-exposed and flexible enough so that an extra backbone ligand from the outer shell can bind to Sr²⁺. In addition, solvent-exposed Ca²⁺ sites with only a few weak charge-donating ligands that can rearrange to fit the strontium's coordination requirements are susceptible to Sr²⁺ displacement. We provide the physical basis of these results and discuss potential novel protein targets of therapeutic Sr²⁺.

General Charge Transfer Dipole Model for AMOEBA-Like Force Fields

The development of highly accurate force fields is always an importance aspect in molecular modeling. In this work, we introduce a general damping-based charge transfer dipole (D-CTD) model to describe the charge transfer energy and the corresponding charge flow for H, C, N, O, P, S, F, Cl, and Br elements in common bio-organic systems. Then, two effective schemes to evaluate the charge flow from the corresponding induced dipole moment between the interacting molecules were also proposed and discussed. The potential applicability of the D-CTD model in ion-containing systems was also demonstrated in a series of ion-water complexes including Li⁺, Na⁺, K⁺, Mg²⁺, Ca²⁺, Fe²⁺, Zn²⁺, Pt²⁺, F^-, Cl^-, Br^-, and I^- ions. In general, the D-CTD model demonstrated good accuracy and good transferability in both charge transfer energy and the corresponding charge flow for a wide range of model systems. By distinguishing the intermolecular charge redistribution (charge transfer) under the influence of an external electric field from the accompanying intramolecular charge redistribution (polarization), the D-CTD model is theoretically consistent with current induced dipole-based polarizable dipole models and hence can be easily implemented and parameterized. Along with our previous work in charge penetration-corrected electrostatics, a bottom-up approach constructed water model was also proposed and demonstrated. The structure-maker and structure-breaker roles of cations and anions were also correctly reproduced using Na⁺, K⁺, Cl^-, and I^- ions in the new water model, respectively. This work demonstrates a cost-effective approach to describe the charge transfer phenomena. The water and ion models also show the feasibility of a modulated development approach for future force fields.

GPDOCK: highly accurate docking strategy for metalloproteins based on geometric probability

Accurately predicting the interaction modes for metalloproteins remains extremely challenging in structure-based drug design and mechanism analysis of enzymatic catalysis due to the complexity of metal coordination in metalloproteins. Here, we report a docking method for metalloproteins based on geometric probability (GPDOCK) with unprecedented accuracy. The docking tests of 10 common metal ions with 9360 metalloprotein-ligand complexes demonstrate that GPDOCK has an accuracy of 94.3% in predicting binding pose. What is more, it can accurately realize the docking of metalloproteins with ligand when one or two water molecules are engaged in the metal ion coordination. Since GPDOCK only depends on the three-dimensional structure of metalloprotein and ligand, structure-based machine learning model is employed for the scoring of binding poses, which significantly improves computational efficiency. The proposed docking strategy can be an effective and efficient tool for drug design and further study of binding mechanism of metalloproteins. The manual of GPDOCK and the code for the logistical regression model used to re-rank the docking results are available at https://github.com/wangkai-zhku/GPDOCK.git.

Perspective on the Current State-of-the-Art of Quantum Computing for Drug Discovery Applications

Computational chemistry is an essential tool in the pharmaceutical industry. Quantum computing is a fast evolving technology that promises to completely shift the computational capabilities in many areas of chemical research by bringing into reach currently impossible calculations. This perspective illustrates the near-future applicability of quantum computation of molecules to pharmaceutical problems. We briefly summarize and compare the scaling properties of state-of-the-art quantum algorithms and provide novel estimates of the quantum computational cost of simulating progressively larger embedding regions of a pharmaceutically relevant covalent protein-drug complex involving the drug Ibrutinib. Carrying out these calculations requires an error-corrected quantum architecture that we describe. Our estimates showcase that recent developments on quantum phase estimation algorithms have dramatically reduced the quantum resources needed to run fully quantum calculations in active spaces of around 50 orbitals and electrons, from estimated over 1000 years using the Trotterization approach to just a few days with sparse qubitization, painting a picture of fast and exciting progress in this nascent field.

From quantum-derived principles underlying cysteine reactivity to combating the COVID-19 pandemic

The COVID‐19 pandemic poses a challenge in coming up with quick and effective means to counter its cause, the SARS‐CoV‐2. Here, we show how the key factors governing cysteine reactivity in proteins derived from combined quantum mechanical/continuum calculations led to a novel multi‐targeting strategy against SARS‐CoV‐2, in contrast to developing potent drugs/vaccines against a single viral target such as the spike protein. Specifically, they led to the discovery of reactive cysteines in evolutionary conserved Zn²⁺‐sites in several SARS‐CoV‐2 proteins that are crucial for viral polypeptide proteolysis as well as viral RNA synthesis, proofreading, and modification. These conserved, reactive cysteines, both free and Zn²⁺‐bound, can be targeted using the same Zn‐ejector drug (disulfiram/ebselen), which enables the use of broad‐spectrum anti‐virals that would otherwise be removed by the virus's proofreading mechanism. Our strategy of targeting multiple, conserved viral proteins that operate at different stages of the virus life cycle using a Zn‐ejector drug combined with other broad‐spectrum anti‐viral drug(s) could enhance the barrier to drug resistance and antiviral effects, as compared to each drug alone. Since these functionally important nonstructural proteins containing reactive cysteines are highly conserved among coronaviruses, our proposed strategy has the potential to tackle future coronaviruses.

This article is categorized under:

Structure and Mechanism > Reaction Mechanisms and Catalysis

Structure and Mechanism > Computational Biochemistry and Biophysics

Electronic Structure Theory > Density Functional Theory

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.

Automatically Constructed Neural Network Potentials for Molecular Dynamics Simulation of Zinc Proteins

The development of accurate and efficient potential energy functions for the molecular dynamics simulation of metalloproteins has long been a great challenge for the theoretical chemistry community. An artificial neural network provides the possibility to develop potential energy functions with both the efficiency of the classical force fields and the accuracy of the quantum chemical methods. In this work, neural network potentials were automatically constructed by using the ESOINN-DP method for typical zinc proteins. For the four most common zinc coordination modes in proteins, the potential energy, atomic forces, and atomic charges predicted by neural network models show great agreement with quantum mechanics calculations and the neural network potential can maintain the coordination geometry correctly. In addition, MD simulation and energy optimization with the neural network potential can be readily used for structural refinement. The neural network potential is not limited by the function form and complex parameterization process, and important quantum effects such as polarization and charge transfer can be accurately considered. The algorithm proposed in this work can also be directly applied to proteins containing other metal ions.

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.

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.

A Fragment Quantum Mechanical Method for Metalloproteins

An accurate energy calculation of metalloprotein is of crucial importance and also a theoretical challenge. In this work, a metal molecular fractionation with conjugate caps (metal-MFCC) approach is developed for efficient linear-scaling quantum calculation of potential energy and atomic forces of metalloprotein. In this approach, the potential energy of a given protein is calculated by a linear combination of potential energies of the neighboring residues, two-body interaction energy between non-neighboring residues that are spatially in close contact and the potential energy of the metal binding group. The calculation of each fragment is embedded in a field of point charges representing the remaining protein environment. Numerical studies were carried out to check the performance of this method, and the calculated potential energies and atomic forces all show excellent agreement with the full system calculations at the M06-2X/6-31G(d) level. By combining the energy calculation with molecular dynamic simulation, we performed an ab initio structural optimization for a zinc finger protein with high efficiency. The present metal-MFCC approach is linear-scaling with a low prefactor and trivially parallelizable. The individual fragment typically contains about 50 atoms, and it is thus possible to be calculated at higher levels of the quantum chemistry method. This fragment method can be routinely applied to perform structural optimization and ab initio molecular dynamic simulation for metalloproteins of any size.

Simulating the Chelate Effect

Despite the rich history of experimental studies focusing on the thermochemistry and kinetics associated with the chelate effect, molecular-level computational studies on the chelate ring opening/ring closure are scarce. The challenge lies in an accurate description of both the metal ion and its aqueous environment. Herein, we demonstrate that an optimized 12-6-4 Lennard-Jones (LJ) model can capture the thermodynamics and provide detailed structural and mechanistic insights into the formation of ethylenediamine (en) complexes with metal ions. The water molecules in the first solvation shell of the metal ion are found to facilitate the chelate ring formation. The optimized parameters further simulate the formation of bis and tris(en) complexes representing the wide applicability of the model to simulate coordination chemistry and self-assembly processes.

The Zinc Linchpin Motif in the DNA Repair Glycosylase MUTYH: Identifying the Zn2+ Ligands and Roles in Damage Recognition and Repair

The DNA base excision repair (BER) glycosylase MUTYH prevents DNA mutations by catalyzing adenine (A) excision from inappropriately formed 8-oxoguanine (8-oxoG):A mismatches. The importance of this mutation suppression activity in tumor suppressor genes is underscored by the association of inherited variants of MUTYH with colorectal polyposis in a hereditary colorectal cancer syndrome known as MUTYH-associated polyposis, or MAP. Many of the MAP variants encompass amino acid changes that occur at positions surrounding the two-metal cofactor-binding sites of MUTYH. One of these cofactors, found in nearly all MUTYH orthologs, is a [4Fe-4S]2+ cluster coordinated by four Cys residues located in the N-terminal catalytic domain. We recently uncovered a second functionally relevant metal cofactor site present only in higher eukaryotic MUTYH orthologs: a Zn2+ ion coordinated by three Cys residues located within the extended interdomain connector (IDC) region of MUTYH that connects the N-terminal adenine excision and C-terminal 8-oxoG recognition domains. In this work, we identified a candidate for the fourth Zn2+ coordinating ligand using a combination of bioinformatics and computational modeling. In addition, using in vitro enzyme activity assays, fluorescence polarization DNA binding assays, circular dichroism spectroscopy, and cell-based rifampicin resistance assays, the functional impact of reduced Zn2+ chelation was evaluated. Taken together, these results illustrate the critical role that the "Zn2+ linchpin motif" plays in MUTYH repair activity by providing for proper engagement of the functional domains on the 8-oxoG:A mismatch required for base excision catalysis. The functional importance of the Zn2+ linchpin also suggests that adjacent MAP variants or exposure to environmental chemicals may compromise Zn2+ coordination, and ability of MUTYH to prevent disease.

Perspective: Quantum mechanical methods in biochemistry and biophysics

In this perspective article, I discuss several research topics relevant to quantum mechanical (QM) methods in biophysical and biochemical applications. Due to the immense complexity of biological problems, the key is to develop methods that are able to strike the proper balance of computational efficiency and accuracy for the problem of interest. Therefore, in addition to the development of novel ab initio and density functional theory based QM methods for the study of reactive events that involve complex motifs such as transition metal clusters in metalloenzymes, it is equally important to develop inexpensive QM methods and advanced classical or quantal force fields to describe different physicochemical properties of biomolecules and their behaviors in complex environments. Maintaining a solid connection of these more approximate methods with rigorous QM methods is essential to their transferability and robustness. Comparison to diverse experimental observables helps validate computational models and mechanistic hypotheses as well as driving further development of computational methodologies.

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.

Balancing the Interactions of Mg2+ in Aqueous Solution and with Nucleic Acid Moieties For a Polarizable Force Field Based on the Classical Drude Oscillator Model

Mg²⁺ ions are important in biological systems, particularly in stabilizing compact RNA folds. Mg²⁺ is strongly polarizing, and representing its interactions in heterogeneous environments is a challenge for empirical force field development. To date, the most commonly used force fields in molecular dynamics simulations utilize a pairwise-additive approximation for electrostatic interactions, which cannot account for the significant polarization response in systems containing Mg²⁺. In the present work, we refine the interactions of Mg²⁺ with water, Cl^- ions, and nucleic acid moieties using a polarizable force field based on the classical Drude oscillator model. By targeting gas-phase quantum mechanical interaction energies and geometries of hydrated complexes, as well as condensed-phase osmotic pressure calculations, we present a model for Mg²⁺ that yields quantitative agreement with experimental measurements of water dissociation free energy and osmotic pressure across a broad range of concentrations. Notable is the direct modeling of steric repulsion between the water Drude oscillators and Mg²⁺ to treat the Pauli exclusion effects associated with overlap of the electron clouds of water molecules in the first hydration shell around Mg²⁺. Combined with the refined interactions with nucleic acid moieties, the present model represents a significant advancement in simulating nucleic acid systems containing Mg²⁺.

Factors Governing the Bridging Water Protonation State in Polynuclear Mg(2+) Proteins

An aqua ligand bridges metal cations in a wide variety of enzymes, many of which are drug targets for various diseases. However, the factors affecting its protonation state and thus biological roles remain elusive. By computing the free energy for replacing the bridging H2O by OH(-) in various model Mg(2+) sites, we have evaluated how the nature of an aqua bridge depends on the site's net charge (i.e., the number of charged ligands in the first and second shell and the number of metal cations), the site's solvent exposure, the ligand's charge-donating ability, the bridging oxygen's hydrogen-bonding interactions, intramolecular proton transfer from the bridging H2O to a nearby carboxylate, and the metal coordination number. The results reveal the key factors dictating the protonation state of bridging H2O and provide guidelines in predicting whether H2O or OH(-) bridges two Mg(2+) in polynuclear sites. This helps to elucidate the nucleophile in the enzyme-catalyzed reaction and the net charge of the metal complex (metal cation and first-shell ligands), which plays a critical role in binding.

Thiol versus hydroxamate as zinc binding group in HDAC inhibition: An ab initio QM/MM molecular dynamics study

Zinc-dependent histone deacetylases (HDACs) play a critical role in transcriptional repression and gene silencing, and are among the most attractive targets for the development of new therapeutics against cancer and various other diseases. Two HDAC inhibitors have been approved by FDA as anti-cancer drugs: one is SAHA whose hydroxamate is directly bound to zinc, the other is FK228 whose active form may use thiol as the zinc binding group. In spite of extensive studies, it remains to be ambiguous regarding how thiol and hydroxamate are bound to the zinc active site of HDACs. In this work, our computational approaches center on Born-Oppenheimer ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics with umbrella sampling, which allow for modeling of the zinc active site with reasonable accuracy while properly including dynamics and effects of protein environment. Meanwhile, an improved short-long effective function (SLEF2) to describe non-bonded interactions between zinc and other atoms has been employed in initial MM equilibrations. Our ab initio QM/MM MD simulations have confirmed that hydroxamate is neutral when it is bound to HDAC8, and found that thiol is deprotonated when directly bound to zinc in the HDAC active site. By comparing thiol and hydroxamate, our results elucidated the differences in their binding environment in the HDAC active sites, and emphasized the importance of the linker design to achieve more specific binding toward class IIa HDACs.

DFTB3 Parametrization for Copper: The Importance of Orbital Angular Momentum Dependence of Hubbard Parameters

We report the parametrization of a density functional tight binding method (DFTB3) for copper in a spin-polarized formulation. The parametrization is consistent with the framework of 3OB for main group elements (ONCHPS) and can be readily used for biological applications that involve copper proteins/peptides. The key to our parametrization is to introduce orbital angular momentum dependence of the Hubbard parameter and its charge derivative, thus allowing the 3d and 4s orbitals to adopt different sizes and responses to the change of charge state. The parametrization has been tested by applying to a fairly broad set of molecules of biological relevance, and the properties of interest include optimized geometries, ligand binding energies, and ligand proton affinities. Compared to the reference QM level (B3LYP/aug-cc-pVTZ, which is shown here to be similar to the B97-1 and CCSD(T) results, in terms of many properties of interest for a set of small copper containing molecules), our parametrization generally gives reliable structural properties for both Cu(I) and Cu(II) compounds, although several exceptions are also noted. For energetics, the results are more accurate for neutral ligands than for charged ligands, likely reflecting the minimal basis limitation of DFTB3; the results generally outperform NDDO based methods such as PM6 and even PBE with the 6-31+G(d,p) basis. For all ligand types, single-point B3LYP calculations at DFTB3 geometries give results very close (∼1-2 kcal/mol) to the reference B3LYP values, highlighting the consistency between DFTB3 and B3LYP structures. Possible further developments of the DFTB3 model for a better treatment of transition-metal ions are also discussed. In the current form, our first generation of DFTB3 copper model is expected to be particularly valuable as a method that drives sampling in systems that feature a dynamical copper binding site.

Systematic Parameterization of Monovalent Ions Employing the Nonbonded Model

Monovalent ions play fundamental roles in many biological processes in organisms. Modeling these ions in molecular simulations continues to be a challenging problem. The 12-6 Lennard-Jones (LJ) nonbonded model is widely used to model monovalent ions in classical molecular dynamics simulations. A lot of parameterization efforts have been reported for these ions with a number of experimental end points. However, some reported parameter sets do not have a good balance between the two Lennard-Jones parameters (the van der Waals (VDW) radius and potential well depth), which affects their transferability. In the present work, via the use of a noble gas curve we fitted in former work (J. Chem. Theory Comput. 2013, 9, 2733), we reoptimized the 12-6 LJ parameters for 15 monovalent ions (11 positive and 4 negative ions) for three extensively used water models (TIP3P, SPC/E, and TIP4P(EW)). Since the 12-6 LJ nonbonded model performs poorly in some instances for these ions, we have also parameterized the 12-6-4 LJ-type nonbonded model (J. Chem. Theory Comput. 2014, 10, 289) using the same three water models. The three derived parameter sets focused on reproducing the hydration free energies (the HFE set) and the ion-oxygen distance (the IOD set) using the 12-6 LJ nonbonded model and the 12-6-4 LJ-type nonbonded model (the 12-6-4 set) overall give improved results. In particular, the final parameter sets showed better agreement with quantum mechanically calculated VDW radii and improved transferability to ion-pair solutions when compared to previous parameter sets.

Ion selectivity strategies of sodium channel selectivity filters

CONSPECTUS: Sodium ion channels selectively transport Na(+) cations across the cell membrane. These integral parts of the cell machinery are implicated in regulating the cardiac, skeletal and smooth muscle contraction, nerve impulses, salt and water homeostasis, as well as pain and taste perception. Their malfunction often results in various channelopathies of the heart, brain, skeletal muscles, and lung; thus, sodium channels are key drug targets for various disorders including cardiac arrhythmias, heart attack, stroke, migraine, epilepsy, pain, cancer, and autoimmune disorders. The ability of sodium channels to discriminate the native Na(+) among other competing ions in the surrounding fluids is crucial for proper cellular functions. The selectivity filter (SF), the narrowest part of the channel's open pore, lined with amino acid residues that specifically interact with the permeating ion, plays a major role in determining Na(+) selectivity. Different sodium channels have different SFs, which vary in the symmetry, number, charge, arrangement, and chemical type of the metal-ligating groups and pore size: epithelial/degenerin/acid-sensing ion channels have generally trimeric SFs lined with three conserved neutral serines and/or backbone carbonyls; eukaryotic sodium channels have EKEE, EEKE, DKEA, and DEKA SFs with an invariant positively charged lysine from the second or third domain; and bacterial voltage-gated sodium (Nav) channels exhibit symmetrical EEEE SFs, reminiscent of eukaryotic voltage-gated calcium channels. How do these different sodium channel SFs achieve high selectivity for Na(+) over its key rivals, K(+) and Ca(2+)? What factors govern the metal competition in these SFs and which of these factors are exploited to achieve Na(+) selectivity in the different sodium channel SFs? The free energies for replacing K(+) or Ca(2+) bound inside different model SFs with Na(+), evaluated by a combination of density functional theory and continuum dielectric calculations, have shed light on these questions. The SFs of epithelial and eukaryotic Nav channels select Na(+) by providing an optimal number and ligating strength of metal ligands as well as a rigid pore whose size fits the cognate Na(+) ideally. On the other hand, the SFs of bacterial Nav channels select Na(+), as the protein matrix attenuates ion-protein interactions relative to ion-solvent interactions by enlarging the pore and allowing water to enter, so the ion interacts indirectly with the conserved glutamates via bridging water molecules. This shows how these various SFs have adapted to the specific physicochemical properties of the native ion, using different strategies to select Na(+) among its contenders.

Parametrization of DFTB3/3OB for magnesium and zinc for chemical and biological applications

We report the parametrization of the approximate density functional theory, DFTB3, for magnesium and zinc for chemical and biological applications. The parametrization strategy follows that established in previous work that parametrized several key main group elements (O, N, C, H, P, and S). This 3OB set of parameters can thus be used to study many chemical and biochemical systems. The parameters are benchmarked using both gas-phase and condensed-phase systems. The gas-phase results are compared to DFT (mostly B3LYP), ab initio (MP2 and G3B3), and PM6, as well as to a previous DFTB parametrization (MIO). The results indicate that DFTB3/3OB is particularly successful at predicting structures, including rather complex dinuclear metalloenzyme active sites, while being semiquantitative (with a typical mean absolute deviation (MAD) of ∼3–5 kcal/mol) for energetics. Single-point calculations with high-level quantum mechanics (QM) methods generally lead to very satisfying (a typical MAD of ∼1 kcal/mol) energetic properties. DFTB3/MM simulations for solution and two enzyme systems also lead to encouraging structural and energetic properties in comparison to available experimental data. The remaining limitations of DFTB3, such as the treatment of interaction between metal ions and highly charged/polarizable ligands, are also discussed.

Parameterization of highly charged metal ions using the 12-6-4 LJ-type nonbonded model in explicit water

Highly charged metal ions act as catalytic centers and structural elements in a broad range of chemical complexes. The nonbonded model for metal ions is extensively used in molecular simulations due to its simple form, computational speed, and transferability. We have proposed and parametrized a 12-6-4 LJ (Lennard-Jones)-type nonbonded model for divalent metal ions in previous work, which showed a marked improvement over the 12-6 LJ nonbonded model. In the present study, by treating the experimental hydration free energies and ion–oxygen distances of the first solvation shell as targets for our parametrization, we evaluated 12-6 LJ parameters for 18 M(III) and 6 M(IV) metal ions for three widely used water models (TIP3P, SPC/E, and TIP4P_EW). As expected, the interaction energy underestimation of the 12-6 LJ nonbonded model increases dramatically for the highly charged metal ions. We then parametrized the 12-6-4 LJ-type nonbonded model for these metal ions with the three water models. The final parameters reproduced the target values with good accuracy, which is consistent with our previous experience using this potential. Finally, tests were performed on a protein system, and the obtained results validate the transferability of these nonbonded model parameters.

Modeling Zn²⁺ release from metallothionein

Mammalian metallothioneins (MTs) comprise a Zn3Cys9 cluster in the β domain and a Zn4Cys11 cluster in the α domain. They play a crucial role in storing and donating Zn(2+) ions to target metalloproteins and have been implicated in several diseases, thus understanding how MTs release Zn(2+) is of widespread interest. In this work, we present a strategy to compute the free energy for releasing Zn(2+) from MTs using a combination of classical molecular dynamics (MD) simulations, quantum-mechanics/molecular-mechanics (QM/MM) minimizations, and continuum dielectric calculations. The methodology is shown to reproduce the experimental observations that (1) the Zn-binding sites do not have equal Zn(2+) affinity and (2) the isolated β domain is thermodynamically less stable and releases Zn(2+) faster with oxidizing agents than the isolated α domain. It was used to compute the free energies for Zn(2+) release from the metal cluster in the absence and presence of the protein matrix (protein architecture and coupled protein-water interactions) to yield the respective disulfide-bonded product. The results show the importance of the protein matrix as well as protein dynamics and coupled conformational changes in accounting for the differential Zn(2+)-releasing propensity of the two domains with oxidizing agents.

Metal-induced conformational changes of human metallothionein-2A: a combined theoretical and experimental study of metal-free and partially metalated intermediates

Electrospray ionization ion mobility mass spectrometry (ESI IM-MS) and molecular dynamics (MD) simulations reveal new insights into metal-induced conformational changes and the mechanism for metalation of human metallothionein-2A (MT), an intrinsically disordered protein. ESI of solutions containing apoMT yields multiple charge states of apoMT; following addition of Cd(2+) to the solution, ESI yields a range of CdiMT (i = 1-7) product ions (see Chen et al. Anal. Chem. 2013, 85, 7826-33). Ion mobility arrival-time distributions (ATDs) for the CdiMT (i = 0-7) ions reveal a diverse population of ion conformations. The ion mobility data clearly show that the conformational diversity for apoMT and partially metalated ions converges toward ordered, compact conformations as the number of bound Cd(2+) ions increase. MD simulations provide additional information on conformation candidates of CdiMT (i = 0-7) that supports the convergence of distinct conformational populations upon metal binding. Integrating the IM-MS and MD data provides a global view that shows stepwise conformational transition of an ensemble as a function of metal ion bound. ApoMT is comprised of a wide range of conformational states that populate between globular-like compact and coil-rich extended conformations. During the initial stepwise metal addition (number of metal ions bound i = 1-3), the metal ions bind to different sites to yield distinct conformations, whereas for i > 4, the conformational changes appear to be domain-specific, attributed to different degrees of disorder of the β domain; the β domain becomes more ordered as additional metal ions are added, promoting convergences to the dumbbell-shaped conformation.

Coordination-directed self-assembly of M12L24 nanocage: effects of kinetic trapping on the assembly process

We demonstrate the spontaneous formation of spherical complex M12L24, which is composed of 12 palladium ions and 24 bidentate ligands, by molecular dynamics simulations. In contrast to our previous study on the smaller M6L8 cage, we found that the larger M12L24 self-assembly process involves noticeable kinetic trapping at lower nuclearity complexes, e.g., M6L12, M8L16, and M9L18. We also found that the kinetic trapping behaviors sensitively depend on the bend angle of ligands and the metal-ligand binding strength. Our results show that these kinetic effects, that have generally been neglected, are important factor in self-assembly structure determination of larger complexes as M12L24 in this study.

Taking into Account the Ion-induced Dipole Interaction in the Nonbonded Model of Ions

Metal ions exist in almost half of the proteins in the protein databank and they serve as structural, electron-transfer and catalytic elements in the metabolic processes of organisms. Molecular Dynamics (MD) simulation is a powerful tool that provides information about biomolecular systems at the atomic level. Coupled with the growth in computing power, algorithms like the Particle Mesh Ewald (PME) method have become the accepted standard when dealing with long-range interactions in MD simulations. The nonbonded model of metal ions consists of an electrostatic plus 12-6 Lennard Jones (LJ) potential and is used largely because of its speed relative to more accurate models. In previous work we found that ideal parameters do not exist that reproduce several experimental properties for M(II) ions simultaneously using the nonbonded model coupled with the PME method due to the underestimation of metal ion-ligand interactions. Via a consideration of the nature of the nonbonded model, we proposed that the observed error largely arises from overlooking charge-induced dipole interactions. The electrostatic plus 12-6 LJ potential model works reasonably well for neutral systems but does struggle with more highly charged systems. In the present work we designed and parameterized a new nonbonded model for metal ions by adding a 1/r4 term to the 12-6 model. We call it the 12-6-4 LJ-type nonbonded model due to its mathematical construction. Parameters were determined for 16 +2 metal ions for the TIP3P, SPC/E and TIP4PEW water models. The final parameters reproduce the experimental hydration free energies (HFE), ion-oxygen distances (IOD) in the first solvation shell and coordination numbers (CN) accurately for the metal ions investigated. Preliminary tests on MgCl2 at different concentrations in aqueous solution and Mg2+--nucleic acid systems show reasonable results suggesting that the present parameters can work in mixed systems. The 12-6-4 LJ-type nonbonded model is readily adopted into standard force fields like AMBER, CHARMM and OPLS-AA with only a modest computational overhead. The new nonbonded model doesn't consider charge-transfer effects explicitly and, hence, may not suitable for the simulation of systems where charge-transfer effects play a decisive role.

A simple but effective modeling strategy for structural properties of non-heme Fe(II) sites in proteins: test of force field models and application to proteins in the AlkB family

To facilitate computational study of proteins in the AlkB family and related α-ketoglutarate/Fe(II)-dependent dioxygenases, we have tested a simple modeling strategy for the non-heme Fe(II) site in which the iron is represented by a simple +2 point charge with Lennard-Jones parameters. Calculations for an AlkB active site model in the gas phase and ∼150 ns molecular dynamics (MD) simulations for two enzyme-dsDNA complexes (E. coli AlkB-dsDNA and ABH2-dsDNA) suggest that this simple modeling strategy provides a satisfactory description of structural properties of the Fe(II) site in AlkB enzymes, provided that care is exercised to control the binding mode of carboxylate (Asp) to the iron. MD simulations using the model for AlkB-dsDNA and ABH2-dsDNA systems find that although the structural features for the latter are overall in good agreement with the crystal structure, the dsDNA, and AlkB-dsDNA interface undergo substantial changes during the MD simulations from the crystal structure. Even for ABH2, new interactions form between a long loop region and dsDNA upon structural relaxation of the loop, supporting the role of this loop in DNA binding despite the lack of interactions between them in the crystal structure. Analysis of DNA backbone torsional distributions helps identify regions that adopt strained conformations. Collectively, the results highlight that crystal packing may have a significant impact on the structure of protein-DNA complexes; the simulations also provide additional insights regarding why AlkB and ABH2 prefer single-strand and double-strand DNA, respectively, as substrate.

A New Quantum Calibrated Force Field for Zinc-Protein Complex

A quantum calibrated polarizable-charge transfer force field (QPCT) has been proposed to accurately describe the interaction dynamics of zinc-protein complexes. The parameters of the QPCT force field were calibrated by quantum chemistry calculation and capture the polarization and charge transfer effect. QPCTs are validated by molecular dynamic simulation of the hydration shell of the zinc ion, five proteins containing the most common zinc-binding sites (ZnCys2His2, ZnCys3His1, ZnCys4, Zn2Cys6), as well as protein-ligand binding energy in zinc protein MMP3. The calculated results show excellent agreement with the experimental measurement and with results from QM/MM simulation, demonstrating that QPCT is accurate enough to maintain the correct structural integrity of the zinc binding pocket and provide accurate interaction dynamics of the zinc-residue complex. The current approach can also be extended to the study of interaction dynamics of other metal-containing proteins by recalibrating the corresponding parameters to the specific complexes.

Structure and dynamics of a dizinc metalloprotein: effect of charge transfer and polarization

Structures and dynamics of a recently designed dizinc metalloprotein (DFsc) (J. Mol. Biol. 2003, 334, 1101) are studied by molecular dynamics simulation using a dynamically adapted polarized force field derived from fragment quantum calculation for protein in solvent. To properly describe the effect of charge transfer and polarization in the present approach, quantum chemistry calculation of the zinc-binding group is periodically performed (on-the-fly) to update the atomic charges of the zinc-binding group during the MD simulation. Comparison of the present result with those obtained from simulations under standard AMBER force field reveals that charge transfer and polarization are critical to maintaining the correct asymmetric metal coordination in the DFsc. Detailed analysis of the result also shows that dynamic fluctuation of the zinc-binding group facilitates solvent interaction with the zinc ions. In particular, the dynamic fluctuation of the zinc-zinc distance is shown to be an important feature of the catalytic function of the di-ion zinc-binding group. Our study demonstrates that the dynamically adapted polarization approach is computationally practical and can be used to study other metalloprotein systems.

A Transferable Non-bonded Pairwise Force Field to Model Zinc Interactions in Metalloproteins

Herein we introduce a novel practical strategy to overcome the well-known challenge of modeling the divalent zinc cation in metalloproteins. The main idea is to design short-long effective functions (SLEF) to describe charge interactions between the zinc ion and all other atoms. This SLEF approach has the following desired features: (1). It is pairwise, additive and compatible with widely used atomic pair-wise force fields for modeling biomolecules; (2). It only changes interactions between the zinc ion and other atoms, and does not affect force field parameters that model other interactions in the system; (3). It is a non-bonded model that is inherently capable to describe different zinc ligands and coordination modes. By optimizing two SLEF parameters as well as zinc vdW parameters through force matching based on Born-Oppenheimer ab initio QM/MM molecular dynamics simulations, we have successfully developed the first SLEF force field (SLEF1) to describe zinc interactions. Extensive molecular dynamics simulations of seven zinc enzyme systems with different coordination ligands and distinct chelation modes (4-,5-,6-fold), including the binuclear zinc active site, yielded zinc coordination numbers and binding distances in good agreement with the corresponding crystal structures as well as ab initio QM/MM MD results. This not only demonstrates the transferability and adequacy of the new SLEF1 force field in describing a variety of zinc proteins, but also indicates that this novel SLEF approach is a promising direction to explore for improving force field description of metal ion interactions.

Polarization effects in molecular mechanical force fields

The focus here is on incorporating electronic polarization into classical molecular mechanical force fields used for macromolecular simulations. First, we briefly examine currently used molecular mechanical force fields and the current status of intermolecular forces as viewed by quantum mechanical approaches. Next, we demonstrate how some components of quantum mechanical energy are effectively incorporated into classical molecular mechanical force fields. Finally, we assess the modeling methods of one such energy component-polarization energy-and present an overview of polarizable force fields and their current applications. Incorporating polarization effects into current force fields paves the way to developing potentially more accurate, though more complex, parameterizations that can be used for more realistic molecular simulations.

Influence of NH-Sgamma bonding interactions on the structure and dynamics of metallothioneins

Mammalian metallothioneins ([Formula: see text]) show a clustered arrangement of the metal ions and a nonregular protein structure. The solution structures of Cd(3)-thiolate cluster containing beta-domain of mouse beta-MT-1 and rat beta-MT-2 show high structural similarities, but widely differing structure dynamics. Molecular dynamics simulations revealed a substantially increased number of NH-Sgamma hydrogen bonds in beta-MT-2, features likely responsible for the increased stability of the Cd(3)-thiolate cluster and the enfolding protein domain. Alterations in the NH-Sgamma hydrogen-bonding network may provide a rationale for the differences in dynamic properties encountered in the beta-domains of MT-1, -2, and -3 isoforms, believed to be essential for their different biological function.