QM/MM X-ray refinement of zinc metalloenzymes

A database overview of metal-coordination distances in metalloproteins

Metalloproteins are ubiquitous in all living organisms and take part in a very wide range of biological processes. For this reason, their experimental characterization is crucial to obtain improved knowledge of their structure and biological functions. The three-dimensional structure represents highly relevant information since it provides insight into the interaction between the metal ion(s) and the protein fold. Such interactions determine the chemical reactivity of the bound metal. The available PDB structures can contain errors due to experimental factors such as poor resolution and radiation damage. A lack of use of distance restraints during the refinement and validation process also impacts the structure quality. Here, the aim was to obtain a thorough overview of the distribution of the distances between metal ions and their donor atoms through the statistical analysis of a data set based on more than 115 000 metal-binding sites in proteins. This analysis not only produced reference data that can be used by experimentalists to support the structure-determination process, for example as refinement restraints, but also resulted in an improved insight into how protein coordination occurs for different metals and the nature of their binding interactions. In particular, the features of carboxylate coordination were inspected, which is the only type of interaction that is commonly present for nearly all metals.

Molecular insights into the binding of carnosine and anserine to human serum carnosinase 1 (CN1)

Carnosine (CAR) and anserine (ANS) are histidine-containing dipeptides that show the therapeutic properties and protective abilities against diabetes and cognitive deficit. Both dipeptides are rich in meat products and have been used as a supplement. However, in humans, both compounds have a short half-life due to the rapid degradation by dizinc carnosinase 1 (CN1) which is a hurdle for its therapeutic application. To date, a comparative study of carnosine- and anserine-CN1 complexes is limited. Thus, in this work, molecular dynamics (MD) simulations were performed to explore the binding of carnosine and anserine to CN1. CN1 comprises 2 chains (Chains A and B). Both monomers are found to work independently and alternatingly. The displacement of Zn2+ pair is found to disrupt the substrate binding. CN1 employs residues from the neighbour chain (H235, T335, and T337) to form the active site. This highlights the importance of a dimer for enzymatic activity. Anserine is more resistant to CN 1 than carnosine because of its bulky and dehydrated imidazole moiety. Although both dipeptides can direct the peptide oxygen to the active Zn2+ which can facilitate the catalytic reaction, the bulky methylated imidazole on anserine promotes various poses that can retard the hydrolytic activity in contrast to carnosine. Anserine is likely to be the temporary competitive inhibitor by retarding the carnosine catabolism.

Small Molecules Promote Selective Denaturation and Degradation of Tubulin Heterodimers through a Low-Barrier Hydrogen Bond

Here, we report a novel mechanism to selectively degrade target proteins. 3-(3-Phenoxybenzyl)amino-β-carboline (PAC), a tubulin inhibitor, promotes selective degradation of αβ-tubulin heterodimers. Biochemical studies have revealed that PAC specifically denatures tubulin, making it prone to aggregation that predisposes it to ubiquitinylation and then degradation. The degradation is mediated by a single hydrogen bond formed between the pyridine nitrogen of PAC and βGlu198, which is identified as a low-barrier hydrogen bond (LBHB). In contrast, another two tubulin inhibitors that only form normal hydrogen bonds with βGlu198 exhibit no degradation effect. Thus, the LBHB accounts for the degradation. We then screened for compounds capable of forming an LBHB with βGlu198 and demonstrated that BML284, a Wnt signaling activator, also promotes tubulin heterodimer degradation through the LBHB. Our study provided a unique example of LBHB function and identified a novel approach to obtain tubulin degraders.

Zinc finger structure determination by NMR: Why zinc fingers can be a handful

Zinc fingers can be loosely defined as protein domains containing one or more tetrahedrally-co-ordinated zinc ions whose role is to stabilise the structure rather than to be involved in enzymatic chemistry; such zinc ions are often referred to as "structural zincs". Although structural zincs can occur in proteins of any size, they assume particular significance for very small protein domains, where they are often essential for maintaining a folded state. Such small structures, that sometimes have only marginal stability, can present particular difficulties in terms of sample preparation, handling and structure determination, and early on they gained a reputation for being resistant to crystallisation. As a result, NMR has played a more prominent role in structural studies of zinc finger proteins than it has for many other types of proteins. This review will present an overview of the particular issues that arise for structure determination of zinc fingers by NMR, and ways in which these may be addressed.

Tris-(2-pyridyl)-pyrazolyl Borate Zinc(II) Complexes: Synthesis, DNA/Protein Binding and In Vitro Cytotoxicity Studies

Zn(II) complexes bearing tris[3-(2-pyridyl)-pyrazolyl] borate (Tppy) ligand (1-3) was synthesized and examined by spectroscopic and analytical tools. Mononuclear [TppyZnCl] (1) has a Zn(II) centre with one arm (pyrazolyl-pyridyl) dangling outside the coordination sphere which is a novel finding in TppyZn(II) chemistry. In complex [TppyZn(H2O)][BF4] (2) hydrogen bonding interaction of aqua moiety stabilizes the dangling arm. In addition, solution state behaviour of complex 1 confirms the tridentate binding mode and reactivity studies show the exogenous axial substituents used to form the [TppyZnN3] (3). The complexes (1-3) were tested for their ability to bind with Calf thymus (CT) DNA and Bovine serum albumin (BSA) wherein they revealed to exhibit good binding constant values with both the biomolecules in the order of 104-105 M-1. The intercalative binding mode with CT DNA was confirmed from the UV-Visible absorption, viscosity, and ethidium bromide (EB) DNA displacement studies. Further, the complexes were tested for in vitro cytotoxic ability on four triple-negative breast cancer (TNBC) cell lines (MDA-MB-231, MDA-MB-468, HCC1937, and Hs 578T). All three complexes (1-3) exhibited good IC50 values (6.81 to 16.87 μM for 24 h as seen from the MTS assay) results which indicated that these complexes were found to be potential anticancer agents against the TNBC cells.

Combining crystallography with quantum mechanics

In standard crystallographic refinement of biomacromolecules, the crystallographic raw data are supplemented by empirical restraints that ensure that the structure makes chemical sense. These restraints are typically accurate for amino acids and nucleic acids, but less so for cofactors, substrates, inhibitors, ligands and metal sites. In quantum refinement, this potential is replaced by more accurate quantum mechanical (QM) calculations. Several implementations have been presented, differing in the level of QM and whether it is used for the entire structure or only for a site of particular interest. It has been shown that the method can improve and correct errors in crystal structures and that it can be used to determine protonation and tautomeric states of various ligands and to decide what is really seen in the structure by refining different interpretations and using standard crystallographic and QM quality measures to decide which fits the structure best.

'All That Glitters Is Not Gold': High-Resolution Crystal Structures of Ligand-Protein Complexes Need Not Always Represent Confident Binding Poses

Our understanding of the structure–function relationships of biomolecules and thereby applying it to drug discovery programs are substantially dependent on the availability of the structural information of ligand–protein complexes. However, the correct interpretation of the electron density of a small molecule bound to a crystal structure of a macromolecule is not trivial. Our analysis involving quality assessment of ~0.28 million small molecule–protein binding site pairs derived from crystal structures corresponding to ~66,000 PDB entries indicates that the majority (65%) of the pairs might need little (54%) or no (11%) attention. Out of the remaining 35% of pairs that need attention, 11% of the pairs (including structures with high/moderate resolution) pose serious concerns. Unfortunately, most users of crystal structures lack the training to evaluate the quality of a crystal structure against its experimental data and, in general, rely on the resolution as a ‘gold standard’ quality metric. Our work aims to sensitize the non-crystallographers that resolution, which is a global quality metric, need not be an accurate indicator of local structural quality. In this article, we demonstrate the use of several freely available tools that quantify local structural quality and are easy to use from a non-crystallographer’s perspective. We further propose a few solutions for consideration by the scientific community to promote quality research in structural biology and applied areas.

Multiscale Quantum Refinement Approaches for Metalloproteins

Biomolecules with metal ion(s) (e.g., metalloproteins) play many important biological roles. However, accurate structural determination of metalloproteins, particularly those containing transition metal ion(s), is challenging due to their complicated electronic structure, complex bonding of metal ions, and high number of conformations in biomolecules. Quantum refinement, which was proposed to combine crystallographic data with computational chemistry methods by several groups, can improve the local structures of some proteins. In this study, a quantum refinement method combining several multiscale computational schemes with experimental (X-ray diffraction) information was developed for metalloproteins. Various quantum refinement approaches using different ONIOM (our own N-layered integrated molecular orbital and molecular mechanics) combinations of quantum mechanics (QM), semiempirical (SE), and molecular mechanics (MM) methods were conducted to assess the performance and reliability on the refined local structure in two metalloproteins. The structures for two (Cu- or Zn-containing) metalloproteins were refined by combining two-layer ONIOM2(QM1/QM2) and ONIOM2(QM/MM) and three-layer ONIOM3(QM1/QM2/MM) schemes with experimental data. The accuracy of the quantum-refined metal binding sites was also examined and compared in these multiscale quantum refinement calculations. ONIOM3(QM/SE/MM) schemes were found to give good results with lower computational costs and were proposed to be a good choice for the multiscale computational scheme for quantum refinement calculations of metal binding site(s) in metalloproteins with high efficiency. Additionally, a two-center ONIOM approach was employed to speed up the quantum refinement calculations for the Zn metalloprotein with two remote active sites/ligands. Moreover, a recent quantum-embedding wavefunction-in-density functional theory (WF-in-DFT) method was also adopted as the high-level method in unprecedented ONIOM2(CCSD-in-B3LYP/MM) and ONIOM3(CCSD-in-B3LYP/SE/MM) calculations, which can be regarded as novel pseudo-three- and pseudo-four-layer ONIOM methods, respectively, to refine the key Zn binding site at the coupled-cluster singles and doubles (CCSD) level. These refined results indicate that multiscale quantum refinement schemes can be used to improve the structural accuracy obtained for local metal binding site(s) in metalloproteins with high efficiency.

Force Field Parameters for Fe2+4S2-4 Clusters of Dihydropyrimidine Dehydrogenase, the 5-Fluorouracil Cancer Drug Deactivation Protein: A Step towards In Silico Pharmacogenomics Studies

The dimeric dihydropyrimidine dehydrogenase (DPD), metalloenzyme, an adjunct anti-cancer drug target, contains highly specialized 4 × Fe²⁺₄S²⁻₄ clusters per chain. These clusters facilitate the catalysis of the rate-limiting step in the pyrimidine degradation pathway through a harmonized electron transfer cascade that triggers a redox catabolic reaction. In the process, the bulk of the administered 5-fluorouracil (5-FU) cancer drug is inactivated, while a small proportion is activated to nucleic acid antimetabolites. The occurrence of missense mutations in DPD protein within the general population, including those of African descent, has adverse toxicity effects due to altered 5-FU metabolism. Thus, deciphering mutation effects on protein structure and function is vital, especially for precision medicine purposes. We previously proposed combining molecular dynamics (MD) and dynamic residue network (DRN) analysis to decipher the molecular mechanisms of missense mutations in other proteins. However, the presence of Fe²⁺₄S²⁻₄ clusters in DPD poses a challenge for such in silico studies. The existing AMBER force field parameters cannot accurately describe the Fe²⁺ center coordination exhibited by this enzyme. Therefore, this study aimed to derive AMBER force field parameters for DPD enzyme Fe²⁺ centers, using the original Seminario method and the collation features Visual Force Field Derivation Toolkit as a supportive approach. All-atom MD simulations were performed to validate the results. Both approaches generated similar force field parameters, which accurately described the human DPD protein Fe²⁺₄S²⁻₄ cluster architecture. This information is crucial and opens new avenues for in silico cancer pharmacogenomics and drug discovery related research on 5-FU drug efficacy and toxicity issues.

QM/MM and molecular dynamics investigation of the mechanism of covalent inhibition of TAK1 kinase

TAK1 is a serine/threonine kinase which is involved in the moderation of cell survival and death via the TNFα signalling pathway. It is also implicated in a range of cancer and anti-inflammatory diseases. Drug discovery efforts on this target have focused on both traditional reversible ATP-binding site inhibitors and increasingly popular irreversible covalent binding inhibitors. Irreversible inhibitors can offer benefits in terms of potency, selectivity and PK/PD meaning they are increasingly pursued where the strategy exists. TAK1 kinase differs from the better-known kinase EGFR in that the reactive cysteine nucleophile targeted by electrophilic inhibitors is located towards the back of the ATP binding site, not at its mouth. While a wealth of structural and computational effort has been spent exploring EGFR, only limited studies on TAK1 have been reported. In this work we report the first QM/MM study on TAK1 aiming to better understand aspects of covalent adduct formation. Our goal is to identify the general base in the catalytic reaction, whether the process proceeds via a stepwise or concerted pathway, and how the highly flexible G-loop and A-loop affect the catalytic cysteine located nearby.

Two similar Schiff-base receptor based quinoline derivate: Highly selective fluorescent probe for Zn(II)

As is known, Zn²⁺ plays a vital role in a variety of biological processes but excessive exposure of Zn²⁺ to human beings can cause toxicity, inducing a series of overt poisoning symptoms and neurodegenerative disorders. Thus, we designed and synthesized two quinoline-derived Schiff-bases HL₁ and HL₂, and investigated the fluorescence emission responses of these two Schiff-bases to various metal ions. A significant enhancement in fluorescence emission band centered at 450 nm was observed in the ethanolic solution of HL₁ with addition of Zn²⁺, while remarkably lower fluorescence emission enhancement was obtained in the case of HL₂ in which one methyl group was introduced to the azomethine carbon. In addition, HL₁ showed good selectivity and high sensitivity towards Zn²⁺ in the existence of other various interfering metal ions, and the reversibility and regeneration of HL₁ were also perfect for extending its applications in environmental and biological systems. Therefore, HL₁ could be identified as a fluorescent probe for sensing Zn²⁺ environmentally and biologically.

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.

User-Friendly Quantum Mechanics: Applications for Drug Discovery

Quantum mechanics (QM) methods provide a fine description of receptor-ligand interactions and of chemical reactions. Their use in drug design and drug discovery is increasing, especially for complex systems including metal ions in the binding sites, for the design of highly selective inhibitors, for the optimization of bi-specific compounds, to understand enzymatic reactions, and for the study of covalent ligands and prodrugs. They are also used for generating molecular descriptors for predictive QSAR/QSPR models and for the parameterization of force fields. Thanks to the continuous increase of computational power offered by GPUs and to the development of sophisticated algorithms, QM methods are becoming part of the standard tools used in computer-aided drug design (CADD). We present the most used QM methods and software packages, and we discuss recent representative applications in drug design and drug discovery.

Theoretical Evaluation of the Reaction Mechanism of Serine Hydroxymethyltransferase

Serine hydroxymethyltransferase (SHMT) is a pyridoxal phosphate (PLP)-dependent enzyme that catalyzes the reversible conversion of serine and tetrahydrofolate (THF) to glycine and 5,10-methylene THF. SHMT is a folate pathway enzyme and is therefore of considerable medical interest due to its role as an important intervention point for antimalarial, anticancer, and antibacterial treatments. Despite considerable experimental effort, the precise reaction mechanism of SHMT remains unclear. In this study, we explore the mechanism of SHMT to determine the roles of active site residues and the nature and the sequence of chemical steps. Molecular dynamics (MD) methods were employed to generate a suitable starting structure which then underwent analysis using hybrid quantum mechanical/molecular mechanical (QM/MM) simulations. The QM region consisted of 12 key residues, two substrates, and explicit solvent. Our results show that the catalytic reaction proceeds according to a retro-aldol synthetic process with His129 acting as the general base in the reaction. The rate-determining step involves the cleavage of the PLP-serine aldimine C_α-C_β bond and the formation of formaldehyde in line with experimental evidence. The pyridyl ring of the PLP-serine aldimine substrate exists in deprotonated form, being stabilized directly by Asp208 via a strong H-bond, as well as through interactions with Arg371, Lys237, and His211, and with the surrounding protein which was electrostatically embedded. This knowledge has the potential to impact the design and development of new inhibitors.

The construction of ZnS–In2S3 nanonests and their heterojunction boosted visible-light photocatalytic/photoelectrocatalytic performance

Nanoscale ZnS–In₂S₃ heterostructure with abundant crystal interfaces delivered excellent visible-light-harvest ability and high photo(electro)catalytic performance.

High-throughput quantum-mechanics/molecular-mechanics (ONIOM) macromolecular crystallographic refinement with PHENIX/DivCon: the impact of mixed Hamiltonian methods on ligand and protein structure

Conventional macromolecular crystallographic refinement relies on often dubious stereochemical restraints, the preparation of which often requires human validation for unusual species, and on rudimentary energy functionals that are devoid of nonbonding effects owing to electrostatics, polarization, charge transfer or even hydrogen bonding. While this approach has served the crystallographic community for decades, as structure-based drug design/discovery (SBDD) has grown in prominence it has become clear that these conventional methods are less rigorous than they need to be in order to produce properly predictive protein-ligand models, and that the human intervention that is required to successfully treat ligands and other unusual chemistries found in SBDD often precludes high-throughput, automated refinement. Recently, plugins to the Python-based Hierarchical ENvironment for Integrated Xtallography (PHENIX) crystallographic platform have been developed to augment conventional methods with the in situ use of quantum mechanics (QM) applied to ligand(s) along with the surrounding active site(s) at each step of refinement [Borbulevych et al. (2014), Acta Cryst D70, 1233-1247]. This method (Region-QM) significantly increases the accuracy of the X-ray refinement process, and this approach is now used, coupled with experimental density, to accurately determine protonation states, binding modes, ring-flip states, water positions and so on. In the present work, this approach is expanded to include a more rigorous treatment of the entire structure, including the ligand(s), the associated active site(s) and the entire protein, using a fully automated, mixed quantum-mechanics/molecular-mechanics (QM/MM) Hamiltonian recently implemented in the DivCon package. This approach was validated through the automatic treatment of a population of 80 protein-ligand structures chosen from the Astex Diverse Set. Across the entire population, this method results in an average 3.5-fold reduction in ligand strain and a 4.5-fold improvement in MolProbity clashscore, as well as improvements in Ramachandran and rotamer outlier analyses. Overall, these results demonstrate that the use of a structure-wide QM/MM Hamiltonian exhibits improvements in the local structural chemistry of the ligand similar to Region-QM refinement but with significant improvements in the overall structure beyond the active site.

Quantum Crystallography: Current Developments and Future Perspectives

Crystallography and quantum mechanics have always been tightly connected because reliable quantum mechanical models are needed to determine crystal structures. Due to this natural synergy, nowadays accurate distributions of electrons in space can be obtained from diffraction and scattering experiments. In the original definition of quantum crystallography (QCr) given by Massa, Karle and Huang, direct extraction of wavefunctions or density matrices from measured intensities of reflections or, conversely, ad hoc quantum mechanical calculations to enhance the accuracy of the crystallographic refinement are implicated. Nevertheless, many other active and emerging research areas involving quantum mechanics and scattering experiments are not covered by the original definition although they enable to observe and explain quantum phenomena as accurately and successfully as the original strategies. Therefore, we give an overview over current research that is related to a broader notion of QCr, and discuss options how QCr can evolve to become a complete and independent domain of natural sciences. The goal of this paper is to initiate discussions around QCr, but not to find a final definition of the field.

Temperature effect on the structure and conformational fluctuations in two zinc knuckles from the mouse mammary tumor virus

Zinc fingers are small protein domains in which zinc plays a structural role, contributing to the stability of the zinc-peptide complex. Zinc fingers are structurally diverse and are present in proteins that perform a broad range of functions in various cellular processes, such as replication and repair, transcription and translation, metabolism and signaling, cell proliferation, and apoptosis. Zinc fingers typically function as interaction modules and bind to a wide variety of compounds, such as nucleic acids, proteins, and small molecules. In this study, we investigated the structural properties, in solution, of the proximal and distal zinc knuckles of the nucleocapsid (NC) protein from the mouse mammary tumor virus (MMTV) (MMTV NC). For this purpose, we performed a series of molecular dynamics simulations in aqueous solution at 300 K, 333 K, and 348 K. The temperature effect was evaluated in terms of root mean square deviation of the backbone atoms and root mean square fluctuation of the coordinating residue atoms. The stability of the zinc coordination sphere was analyzed based upon the time profile of the interatomic distances between the zinc ions and the chelator atoms. The results indicate that the hydrophobic character of the proximal zinc finger is dominant at 333 K. The low mobility of the coordinating residues suggests that the strong electrostatic effect exerted by the zinc ion on its coordinating residues is not influenced by the increase in temperature. The evolution of the structural parameters of the coordination sphere of the distal zinc finger at 300 K gives us a reasonable picture of the unfolding pathway, as proposed by Bombarda and coworkers (Bombarda et al., 2005), which can predict the binding order of the four conserved ligand-binding residues. Our results support the conclusion that the structural features can vary significantly between the two zinc knuckles of MMTV NC.

Quantum crystallography

Approximate wavefunctions can be improved by constraining them to reproduce observations derived from diffraction and scattering experiments. Conversely, charge density models, incorporating electron-density distributions, atomic positions and atomic motion, can be improved by supplementing diffraction experiments with quantum chemically calculated, tailor-made electron densities (form factors). In both cases quantum chemistry and diffraction/scattering experiments are combined into a single, integrated tool. The development of quantum crystallographic research is reviewed. Some results obtained by quantum crystallography illustrate the potential and limitations of this field.

QM/MM investigation of the reaction rates of substrates of 2,3-dimethylmalate lyase: A catabolic protein isolated from Aspergillus niger

Aspergillus niger is an industrially important microorganism used in the production of citric acid. It is a common cause of food spoilage and represents a health issue for patients with compromised immune systems. Recent studies on Aspergillus niger have revealed details on the isocitrate lyase (ICL) superfamily and its role in catabolism, including (2R, 3S)-dimethylmalate lyase (DMML). Members of this and related lyase super families are of considerable interest as potential treatments for bacterial and fungal infections, including Tuberculosis. In our efforts to better understand this class of protein, we investigate the catalytic mechanism of DMML, studying five different substrates and two different active site metals configurations using molecular dynamics (MD) and hybrid quantum mechanics/molecular mechanics (QM/MM) calculations. We show that the predicted barriers to reaction for the substrates show good agreement with the experimental kcat values. This results help to confirm the validity of the proposed mechanism and open up the possibility of developing novel mechanism based inhibitors specifically for this target.

QM/MM investigation of the catalytic mechanism of angiotensin-converting enzyme

Angiotensin-converting enzyme (ACE) converts angiotensin I to angiotensin II and degrades bradykinin and other vasoactive peptides. ACE inhibitors are used to treat diseases such as hypertension and heart failure. It is thus highly desirable to understand the catalytic mechanism of ACE, as this should facilitate the design of more powerful and selective ACE inhibitors. ACE exhibits two different active domains, the C-domain and the N-domain. In this work, we systematically investigated the inhibitor- and substrate-binding patterns in the N-domain of human ACE using a combined quantum mechanical and molecular mechanical approach. The hydrolysis of hippuryl-histidyl-leucine (HHL) as catalyzed by the N-domain of human somatic ACE was explored, and the effects of chloride ion on the overall reaction were also investigated. Two models, one with and one without a chloride ion at the first binding position, were then designed to examine the chloride dependence of inhibitor-substrate binding and the catalytic mechanism. Our calculations indicate that the hydrolysis reaction follows a stepwise general base/general acid catalysis path. The estimated mean free energy barrier height in the two models is about 15.6 kcal/mol, which agrees very well with the experimentally estimated value of 15.8 kcal/mol. Our simulations thus suggest that the N-domain is in a mixed form during ACE-catalyzed hydrolysis, with the single-chloride-ion and the double-chloride-ion forms existing simultaneously. Graphical Abstract Superposition of ACE C- and N- domains.

Exploration of the zinc finger motif in controlling activity of matrix metalloproteinases

Discovering ways to control the activity of matrix metalloproteinases (MMPs), zinc-dependent enzymes capable of degrading extracellular matrix proteins, is an important field of cancer research. We report here a novel strategy for assembling MMP inhibitors on the basis of oligopeptide ligands by exploring the pattern known as the zinc finger motif. Advanced molecular modeling tools were used to characterize the structural binding motifs of experimentally tested MMP inhibitors, as well as those of newly proposed peptidomimetics, in their zinc-containing active sites. The results of simulations based on the quantum mechanics/molecular mechanics (QM/MM) approach and Car–Parrinello molecular dynamics with QM/MM potentials demonstrate that, upon binding of Regasepin1, a known MMP-9 inhibitor, the Zn²⁺(His₃) structural element is rearranged to the Zn²⁺(Cys₂His₂) zinc finger motif, in which two Cys residues are borrowed from the ligand. Following consideration of the crystal structure of MMP-2 with its inhibitor, the oligopeptide APP-IP, we proposed a new peptidomimetic with two replacements in the substrate, Tyr3Cys and Asp6Cys. Simulations show that this peptide variant blocks an enzyme active site by the Zn²⁺(Cys₂His₂) zinc finger construct. Similarly, a natural substrate of MMP-2, Ace-Gln-Gly ∼ Ile-Ala-Gly-Nme, can be converted to an inhibiting compound by two replacements, Ile by Cys and Gly by the d isomer of Cys, favoring formation of the zinc finger motif.

The dynamics of zinc sites in proteins: electronic basis for coordination sphere expansion at structural sites

The functional role assumed by zinc in proteins is closely tied to the variable dynamics around its coordination sphere arising by virtue of its flexibility in bonding. Modern experimental and computational methods allow the detection and study of previously unknown features of bonding between zinc and its ligands in protein environment. These discoveries are occurring just in time as novel biological functions of zinc, which involve rather unconventional coordination trends, are emerging. In this sense coordination sphere expansion of structural zinc sites, as observed in our previous experiments, is a novel phenomenon. Here we explore the electronic and structural requirements by simulating this phenomenon in structural zinc sites using DFT computations. For this purpose, we have chosen MPW1PW91 and a mixed basis set combination as the DFT method through benchmarking, because it accurately reproduces structural parameters of experimentally characterized zinc compounds. Using appropriate models, we show that the greater ionic character of zinc coordination would allow for coordination sphere expansion if the steric and electrostatic repulsions of the ligands are attenuated properly. Importantly, through the study of electronic and structural aspects of the models used, we arrive at a comprehensive bonding model, explaining the factors that influence coordination of zinc in proteins. The proposed model along with the existing knowledge would enhance our ability to predict zinc binding sites in proteins, which is today of growing importance given the predicted enormity of the zinc proteome.

Using quantum mechanical approaches to study biological systems

Conspectus Quantum mechanics (QM) has revolutionized our understanding of the structure and reactivity of small molecular systems. Given the tremendous impact of QM in this research area, it is attractive to believe that this could also be brought into the biological realm where systems of a few thousand atoms and beyond are routine. Applying QM methods to biological problems brings an improved representation to these systems by the direct inclusion of inherently QM effects such as polarization and charge transfer. Because of the improved representation, novel insights can be gleaned from the application of QM tools to biomacromolecules in aqueous solution. To achieve this goal, the computational bottlenecks of QM methods had to be addressed. In semiempirical theory, matrix diagonalization is rate limiting, while in density functional theory or Hartree-Fock theory electron repulsion integral computation is rate-limiting. In this Account, we primarily focus on semiempirical models where the divide and conquer (D&C) approach linearizes the matrix diagonalization step with respect to the system size. Through the D&C approach, a number of applications to biological problems became tractable. Herein, we provide examples of QM studies on biological systems that focus on protein solvation as viewed by QM, QM enabled structure-based drug design, and NMR and X-ray biological structure refinement using QM derived restraints. Through the examples chosen, we show the power of QM to provide novel insights into biological systems, while also impacting practical applications such as structure refinement. While these methods can be more expensive than classical approaches, they make up for this deficiency by the more realistic modeling of the electronic nature of biological systems and in their ability to be broadly applied. Of the tools and applications discussed in this Account, X-ray structure refinement using QM models is now generally available to the community in the refinement package Phenix. While the power of this approach is manifest, challenges still remain. In particular, QM models are generally applied to static structures, so ways in which to include sampling is an ongoing challenge. Car-Parrinello or Born-Oppenheimer molecular dynamics approaches address the short time scale sampling issue, but how to effectively use QM to study phenomenon covering longer time scales will be the focus of future research. Finally, how to accurately and efficiently include electron correlation effects to facilitate the modeling of, for example, dispersive interactions, is also a major hurdle that a broad range of groups are addressing The use of QM models in biology is in its infancy, leading to the expectation that the most significant use of these tools to address biological problems will be seen in the coming years. It is hoped that while this Account summarizes where we have been, it will also help set the stage for future research directions at the interface of quantum mechanics and biology.

Accurate macromolecular crystallographic refinement: incorporation of the linear scaling, semiempirical quantum-mechanics program DivCon into the PHENIX refinement package

Macromolecular crystallographic refinement relies on sometimes dubious stereochemical restraints and rudimentary energy functionals to ensure the correct geometry of the model of the macromolecule and any covalently bound ligand(s). The ligand stereochemical restraint file (CIF) requires a priori understanding of the ligand geometry within the active site, and creation of the CIF is often an error-prone process owing to the great variety of potential ligand chemistry and structure. Stereochemical restraints have been replaced with more robust functionals through the integration of the linear-scaling, semiempirical quantum-mechanics (SE-QM) program DivCon with the PHENIX X-ray refinement engine. The PHENIX/DivCon package has been thoroughly validated on a population of 50 protein-ligand Protein Data Bank (PDB) structures with a range of resolutions and chemistry. The PDB structures used for the validation were originally refined utilizing various refinement packages and were published within the past five years. PHENIX/DivCon does not utilize CIF(s), link restraints and other parameters for refinement and hence it does not make as many a priori assumptions about the model. Across the entire population, the method results in reasonable ligand geometries and low ligand strains, even when the original refinement exhibited difficulties, indicating that PHENIX/DivCon is applicable to both single-structure and high-throughput crystallography.

Force field independent metal parameters using a nonbonded dummy model

The cationic dummy atom approach provides a powerful nonbonded description for a range of alkaline-earth and transition-metal centers, capturing both structural and electrostatic effects. In this work we refine existing literature parameters for octahedrally coordinated Mn²⁺, Zn²⁺, Mg²⁺, and Ca²⁺, as well as providing new parameters for Ni²⁺, Co²⁺, and Fe²⁺. In all the cases, we are able to reproduce both M²⁺–O distances and experimental solvation free energies, which has not been achieved to date for transition metals using any other model. The parameters have also been tested using two different water models and show consistent performance. Therefore, our parameters are easily transferable to any force field that describes nonbonded interactions using Coulomb and Lennard-Jones potentials. Finally, we demonstrate the stability of our parameters in both the human and Escherichia coli variants of the enzyme glyoxalase I as showcase systems, as both enzymes are active with a range of transition metals. The parameters presented in this work provide a valuable resource for the molecular simulation community, as they extend the range of metal ions that can be studied using classical approaches, while also providing a starting point for subsequent parametrization of new metal centers.

Computational characterization of ketone-ketal transformations at the active site of matrix metalloproteinases

We modeled the first steps of hydrolysis reactions of a natural oligopeptide substrate of matrix metalloproteinase MMP-2 as well as of a substrate analogue. In the latter, the scissile amide group is substituted by a ketomethylene group which can be transformed to the ketal group upon binding of this compound to the enzyme active site. According to our quantum mechanical-molecular mechanical (QM/MM) calculations, the reaction of the ketone-ketal transformation proceeds with a low energy barrier (3.4 kcal/mol) and a high equilibrium constant (10(4)). The reaction product with the ketal group formed directly at the active site of the enzyme works as an inhibitor that chelates the zinc ion. On the other hand, the oligopeptide mimetic retains molecular groups responsible for binding of this compound to the enzyme active site. This example illustrates a strategy to design MMP inhibitors in situ by using data on binding specificity of substrates to a particular type of MMP and details of the reaction mechanism.

Conformational analysis and parallel QM/MM X-ray refinement of protein bound anti-Alzheimer drug donepezil

The recognition and association of donepezil with acetylcholinesterase (AChE) has been extensively studied in the past several decades because of the former's use as a palliative treatment for mild Alzheimer disease. Herein we examine the conformational properties of donepezil and we re-examine the donepezil-AChE crystal structure using combined quantum mechanical/molecular mechanical (QM/MM) X-ray refinement tools. Donepezil's conformational energy surface was explored using the M06 suite of density functionals and with the MP2/complete basis set (CBS) method using the aug-cc-pVXZ (X = D and T) basis sets. The donepezil-AChE complex (PDB 1EVE) was also re-refined through a parallel QM/MM X-ray refinement approach based on an in-house ab initio code QUICK, which uses the message passing interface (MPI) in a distributed SCF algorithm to accelerate the calculation via parallelization. In the QM/MM re-refined donepezil structure, coordinate errors that previously existed in the PDB deposited geometry were improved leading to an improvement of the modeling of the interaction between donepezil and the aromatic side chains located in the AChE active site gorge. As a result of the re-refinement there was a 93% reduction in the donepezil conformational strain energy versus the original PDB structure. The results of the present effort offer further detailed structural and biochemical inhibitor-AChE information for the continued development of more effective and palliative treatments of Alzheimer disease.

Solution NMR refinement of a metal ion bound protein using metal ion inclusive restrained molecular dynamics methods

Correctly calculating the structure of metal coordination sites in a protein during the process of nuclear magnetic resonance (NMR) structure determination and refinement continues to be a challenging task. In this study, we present an accurate and convenient means by which to include metal ions in the NMR structure determination process using molecular dynamics (MD) simulations constrained by NMR-derived data to obtain a realistic and physically viable description of the metal binding site(s). This method provides the framework to accurately portray the metal ions and its binding residues in a pseudo-bond or dummy-cation like approach, and is validated by quantum mechanical/molecular mechanical (QM/MM) MD calculations constrained by NMR-derived data. To illustrate this approach, we refine the zinc coordination complex structure of the zinc sensing transcriptional repressor protein Staphylococcus aureus CzrA, generating over 130 ns of MD and QM/MM MD NMR-data compliant sampling. In addition to refining the first coordination shell structure of the Zn(II) ion, this protocol benefits from being performed in a periodically replicated solvation environment including long-range electrostatics. We determine that unrestrained (not based on NMR data) MD simulations correlated to the NMR data in a time-averaged ensemble. The accurate solution structure ensemble of the metal-bound protein accurately describes the role of conformational sampling in allosteric regulation of DNA binding by zinc and serves to validate our previous unrestrained MD simulations of CzrA. This methodology has potentially broad applicability in the structure determination of metal ion bound proteins, protein folding and metal template protein-design studies.

Combined quantum mechanics/molecular mechanics (QM/MM) methods in computational enzymology

Computational enzymology is a rapidly maturing field that is increasingly integral to understanding mechanisms of enzyme-catalyzed reactions and their practical applications. Combined quantum mechanics/molecular mechanics (QM/MM) methods are important in this field. By treating the reacting species with a quantum mechanical method (i.e., a method that calculates the electronic structure of the active site) and including the enzyme environment with simpler molecular mechanical methods, enzyme reactions can be modeled. Here, we review QM/MM methods and their application to enzyme-catalyzed reactions to investigate fundamental and practical problems in enzymology. A range of QM/MM methods is available, from cheaper and more approximate methods, which can be used for molecular dynamics simulations, to highly accurate electronic structure methods. We discuss how modeling of reactions using such methods can provide detailed insight into enzyme mechanisms and illustrate this by reviewing some recent applications. We outline some practical considerations for such simulations. Further, we highlight applications that show how QM/MM methods can contribute to the practical development and application of enzymology, e.g., in the interpretation and prediction of the effects of mutagenesis and in drug and catalyst design.

Evaluating the enthalpic contribution to ligand binding using QM calculations: effect of methodology on geometries and interaction energies

As a result of research on ligand efficiency in the pharmaceutical industry, there is greater focus on optimizing the strength of polar interactions within receptors, so that the contribution of overall size and lipophilicity to binding can be decreased. A number of quantum mechanical (QM) methods involving simple probes are available to assess the H-bonding potential of different heterocycles or functional groups. However, in most receptors, multiple features are present, and these have distinct directionality, meaning very minimalist models may not be so ideal to describe the interactions. We describe how the use of gas phase QM models of kinase protein-ligand complex, which can more closely mimic the polar features of the active site region, can prove useful in assessing alterations to a core template, or different substituents. We investigate some practical issues surrounding the use of QM cluster models in structure based design (SBD). These include the choice of the method; semi-empirical, density functional theory or ab-initio, the choice of the basis set, whether to include implicit or explicit solvation, whether BSSE should be included, etc. We find a combination of the M06-2X method and the 6-31G* basis set is sufficiently rapid, and accurate, for the computation of structural and energetic parameters for this system.

The central role of Gln63 for the hydrogen bonding network and UV-visible spectrum of the AppA BLUF domain

In blue-light sensing using flavin (BLUF) domains, the side-chain orientation of key residues close to the flavin chromophore is still under debate. We report quantum refinements of the wild-type AppA BLUF protein from Rhodobacter sphaeroides starting from two published X-ray structures (1YRX and 2IYG) with different arrangements of the residues around the chromophore. Quantum refinement uses the same experimental X-ray raw data as conventional refinement, but includes data from quantum mechanics/molecular mechanics (QM/MM) calculations as restraints, which is expected to be more reliable than the normally employed MM data. In addition to quantum refinement, pure QM/MM geometry optimizations are performed for the 1YRX and 2IYG structures and for five models derived therefrom. Vertical excitation energies are computed at the QM(DFT/MRCI)/MM level to assess the resulting structures. The experimental absorption maximum of the dark state of wild-type AppA is well reproduced for structures that contain the Gln63 residue in 1YRX-type orientation. The computed excitation energies are red-shifted for structures with a flipped Gln63 residue in 2IYG-type orientation. The calculated 1YRX- and 2IYG-type hydrogen-bonding networks are discussed in detail, particularly with regard to the orientation of the chromophore and the Gln63, Trp104, and Met106 residues.

Conformational Analysis of Free and Bound Retinoic Acid

The conformational profiles of unbound all-trans and 9-cis retinoic acid (RA) have been determined using classical and quantum mechanical calculations. Sixty-six all-trans-RA (ATRA) and forty-eight 9-cis-RA energy minimum conformers were identified via HF/6-31G* geometry optimizations in vacuo. Their relative conformational energies were estimated utilizing the M06, M06-2x and MP2 methods combined with the 6-311+G(d,p), aug-cc-pVDZ and aug-cc-pVTZ basis sets, as well as complete basis set MP2 extrapolations using the latter two basis sets. Single-point energy calculations performed with the M06-2x density functional were found to yield similar results to MP2/CBS for the low-energy retinoic acid conformations. Not unexpectedly, the conformational propensities of retinoic acid were governed by the orientation and arrangement of the torsion angles associated with the polyene tail. We also used previously reported QM/MM X-ray refinement results on four ATRA-protein crystal structures plus one newly refined 9-cis-RA complex (PDB ID 1XDK) in order to investigate the conformational preferences of bound retinoic acid. In the re-refined RA conformers the conjugated double bonds are nearly coplanar, which is consistent with the global minimum identified by the Omega/QM method rather than the corresponding crystallographically determined conformations given in the PDB. Consequently, a 91.3% average reduction of the local strain energy in the gas phase, as well as 92.1% in PCM solvent, was observed using the QM/MM refined structures versus the PDB deposited RA conformations. These results thus demonstrate that our QM/MM X-ray refinement approach can significantly enhance the quality of X-ray crystal structures refined by conventional refinement protocols, thereby providing reliable drug-target structural information for use in structure-based drug discovery applications.

PDB ligand conformational energies calculated quantum-mechanically

We present here a greatly updated version of an earlier study on the conformational energies of protein-ligand complexes in the Protein Data Bank (PDB) [Nicklaus et al. Bioorg. Med. Chem. 1995, 3, 411-428], with the goal of improving on all possible aspects such as number and selection of ligand instances, energy calculations performed, and additional analyses conducted. Starting from about 357,000 ligand instances deposited in the 2008 version of the Ligand Expo database of the experimental 3D coordinates of all small-molecule instances in the PDB, we created a "high-quality" subset of ligand instances by various filtering steps including application of crystallographic quality criteria and structural unambiguousness. Submission of 640 Gaussian 03 jobs yielded a set of about 415 successfully concluded runs. We used a stepwise optimization of internal degrees of freedom at the DFT level of theory with the B3LYP/6-31G(d) basis set and a single-point energy calculation at B3LYP/6-311++G(3df,2p) after each round of (partial) optimization to separate energy changes due to bond length stretches vs bond angle changes vs torsion changes. Even for the most "conservative" choice of all the possible conformational energies-the energy difference between the conformation in which all internal degrees of freedom except torsions have been optimized and the fully optimized conformer-significant energy values were found. The range of 0 to ~25 kcal/mol was populated quite evenly and independently of the crystallographic resolution. A smaller number of "outliers" of yet higher energies were seen only at resolutions above 1.3 Å. The energies showed some correlation with molecular size and flexibility but not with crystallographic quality metrics such as the Cruickshank diffraction-component precision index (DPI) and R(free)-R, or with the ligand instance-specific metrics such as occupancy-weighted B-factor (OWAB), real-space R factor (RSR), and real-space correlation coefficient (RSCC). We repeated these calculations with the solvent model IEFPCM, which yielded energy differences that were generally somewhat lower than the corresponding vacuum results but did not produce a qualitatively different picture. Torsional sampling around the crystal conformation at the molecular mechanics level using the MMFF94s force field typically led to an increase in energy.

Structure modeling of a metalloendopeptidase from Corynebacterium pseudotuberculosis

Metalloendopeptidases are zinc-dependent hydrolases enzymes with many different roles in biological systems, ranging from remodeling conjunctive tissue to removing signaling sequences from nascent proteins. Here, we describe the three-dimensional structure of the metalloendopeptidase from Corynebacterium pseudotuberculosis generated by homology modeling and molecular dynamics. Analysis of key distances shows that His-132, Asp-136, His-211, Leu-212 and one molecule of water play an important role in the protein-Zn(2+) ion interaction. The model obtained may provide structural insights into this enzyme and can be useful for the design of new caseous lymphadenitis vaccines based on genetic attenuation from key point mutation.

QM/MM refinement and analysis of protein bound retinoic acid

Retinoic acid (RA) is a vitamin A derivative, which modifies the appearance of fine wrinkles and roughness of facial skin and treats acne and activates gene transcription by binding to heterodimers of the retinoic acid receptor (RAR) and the retinoic X receptor (RXR). There are series of protein bound RA complexes available in the protein databank (PDB), which provides a broad range of information about the different bioactive conformations of RA. To gain further insights into the observed bioactive RA conformations we applied quantum mechanic (QM)/molecular mechanic (MM) approaches to re-refine the available RA protein-ligand complexes. MP2 complete basis set (CBS) extrapolations single energy calculations are also carried out for both the experimental conformations and QM optimized geometries of RA in the gas as well as solution phase. The results demonstrate that the re-refined structures show better geometries for RA than seen in the originally deposited PDB structures through the use of QMs for the ligand in the X-ray refinement procedure. QM/MM re-refined conformations also reduced the computed strain energies found in the deposited crystal conformations for RA. Finally, the dependence of ligand strain on resolution is analyzed. It is shown that ligand strain is not converged in our calculations and is likely an artifact of the typical resolutions employed to study protein-ligand complexes.

CSAR benchmark exercise of 2010: combined evaluation across all submitted scoring functions

As part of the Community Structure-Activity Resource (CSAR) center, a set of 343 high-quality, protein–ligand crystal structures were assembled with experimentally determined K_d or K_i information from the literature. We encouraged the community to score the crystallographic poses of the complexes by any method of their choice. The goal of the exercise was to (1) evaluate the current ability of the field to predict activity from structure and (2) investigate the properties of the complexes and methods that appear to hinder scoring. A total of 19 different methods were submitted with numerous parameter variations for a total of 64 sets of scores from 16 participating groups. Linear regression and nonparametric tests were used to correlate scores to the experimental values. Correlation to experiment for the various methods ranged R² = 0.58–0.12, Spearman ρ = 0.74–0.37, Kendall τ = 0.55–0.25, and median unsigned error = 1.00–1.68 pK_d units. All types of scoring functions—force field based, knowledge based, and empirical—had examples with high and low correlation, showing no bias/advantage for any particular approach. The data across all the participants were combined to identify 63 complexes that were poorly scored across the majority of the scoring methods and 123 complexes that were scored well across the majority. The two sets were compared using a Wilcoxon rank-sum test to assess any significant difference in the distributions of >400 physicochemical properties of the ligands and the proteins. Poorly scored complexes were found to have ligands that were the same size as those in well-scored complexes, but hydrogen bonding and torsional strain were significantly different. These comparisons point to a need for CSAR to develop data sets of congeneric series with a range of hydrogen-bonding and hydrophobic characteristics and a range of rotatable bonds.

Inhibitor and substrate binding by angiotensin-converting enzyme: quantum mechanical/molecular mechanical molecular dynamics studies

Angiotensin-converting enzyme (ACE) is an important zinc-dependent hydrolase responsible for converting the inactive angiotensin I to the vasoconstrictor angiotensin II and for inactivating the vasodilator bradykinin. However, the substrate binding mode of ACE has not been completely understood. In this work, we propose a model for an ACE Michaelis complex based on two known X-ray structures of inhibitor-enzyme complexes. Specifically, the human testis angiotensin-converting enzyme (tACE) complexed with two clinic drugs were first investigated using a combined quantum mechanical and molecular mechanical (QM/MM) approach. The structural parameters obtained from the 550 ps molecular dynamics simulations are in excellent agreement with the X-ray structures, validating the QM/MM approach. Based on these structures, a model for the Michaelis complex was proposed and simulated using the same computational protocol. Implications to ACE catalysis are discussed.

Polarizable Atomic Multipole X-Ray Refinement: Particle Mesh Ewald Electrostatics for Macromolecular Crystals

Refinement of macromolecular models from X-ray crystallography experiments benefits from prior chemical knowledge at all resolutions. As the quality of the prior chemical knowledge from quantum or classical molecular physics improves, in principle so will resulting structural models. Due to limitations in computer performance and electrostatic algorithms, commonly used macromolecules X-ray crystallography refinement protocols have had limited support for rigorous molecular physics in the past. For example, electrostatics is often neglected in favor of nonbonded interactions based on a purely repulsive van der Waals potential. In this work we present advanced algorithms for desktop workstations that open the door to X-ray refinement of even the most challenging macromolecular data sets using state-of-the-art classical molecular physics. First we describe theory for particle mesh Ewald (PME) summation that consistently handles the symmetry of all 230 space groups, replicates of the unit cell such that the minimum image convention can be used with a real space cutoff of any size and the combination of space group symmetry with replicates. An implementation of symmetry accelerated PME for the polarizable atomic multipole optimized energetics for biomolecular applications (AMOEBA) force field is presented. Relative to a single CPU core performing calculations on a P1 unit cell, our AMOEBA engine called Force Field X (FFX) accelerates energy evaluations by more than a factor of 24 on an 8-core workstation with a Tesla GPU coprocessor for 30 structures that contain 240 000 atoms on average in the unit cell. The benefit of AMOEBA electrostatics evaluated with PME for macromolecular X-ray crystallography refinement is demonstrated via rerefinement of 10 crystallographic data sets that range in resolution from 1.7 to 4.5 Å. Beginning from structures obtained by local optimization without electrostatics, further optimization using AMOEBA with PME electrostatics improved agreement of the model with the data (Rfree was lowered by 0.5%), improved geometric features such as favorable (ϕ, ψ) backbone conformations, and lowered the average potential energy per residue by over 10 kcal/mol. Furthermore, the MolProbity structure validation tool indicates that the geometry of these rerefined structures is consistent with X-ray crystallographic data collected up to 2.2 Å, which is 0.9 Å better than the actual mean quality (3.1 Å). We conclude that polarizable AMOEBA-assisted X-ray refinement offers advantages to methods that neglect electrostatics and is now efficient enough for routine use.

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.