Molecular dynamics simulations of metalloproteins

Applications of the Newly Developed Force-Field Parameters Uncover a Dynamic Nature of Ω-Loop C in the Lys-Ligated Alkaline Form of Cytochrome c

Lys-ligated cytochromes make up an emerging family of heme proteins. Density functional theory calculations on the amine/imidazole-ligated c-type ferric heme were employed to develop force-field parameters for molecular dynamics (MD) simulations of structural and dynamic features of these proteins. The new force-field parameters were applied to the alkaline form of yeast iso-1 cytochrome c to rationalize discrepancies resulting from distinct experimental conditions in prior structural studies and to provide insights into the mechanisms of the alkaline transition. Our simulations have revealed the dynamic nature of Ω-loop C in the Lys-ligated protein and its unfolding in the Lys-ligated conformer having this loop in the same position as in the native Met-ligated protein. The proximity of Tyr67 or Tyr74 to the Lys ligand of ferric heme iron suggests a possible mechanism of the backward alkaline transition where a proton donor Tyr assists in Lys dissociation. The developed force-field parameters will be useful in structural and dynamic characterization of other native or engineered Lys-ligated heme proteins.

Origin of Enantioselectivity in Engineered Cytochrome c-Catalyzed Carbon-Radical FePP Hydrolysis Revealed Using QM/MM (ABEEM Polarizable Force Field) and MD Simulations

The origin of highly efficient asymmetric aminohydroxylation of styrene catalyzed by engineered cytochrome c is investigated by the developed Atom-Bond Electronegativity Equalization Method polarizable force field (ABEEM PFF), which is a combined outcome of electronic and steric effects. Model molecules were used to establish the charge parameters of the ABEEM PFF, for which the bond-stretching and angle-bending parameters were obtained by using a combination of modified Seminario and scan methods. The interactions between carbon-radical Fe-porphyrin (FePP) and waters are simulated by molecular dynamics, which shows a clear preference for the pre-R over the pre-S. This preference is attributed to the hydrogen-bond between the mutated 100S and 101P residues as well as van der Waals interactions, enforcing a specific conformation of the carbon-radical FePP complex within the binding pocket. Meanwhile, the hydrogen-bond between water and the nitrogen atom in the active intermediate dictates the stereochemical outcome. Quantum mechanics/molecular mechanics (QM/MM (ABEEM PFF)) and free-energy perturbation calculations elucidate that the 3RTS is characterized by sandwich-like structure among adjacent amino acid residues, which exhibits greater stability than crowed arrangement in 3STS and enables the R enantiomer to form more favorably. Thus, this study provides mechanistic insight into the catalytic reaction of hemoproteins.

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

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

Role of Metal Cofactor in Enhanced Thermal Stability of Azurin

Metal cofactors are critical centers for different biochemical processes of metalloproteins, and often, this metal coordination renders additional structural stability. In this study, we explore the additional stability conferred by the copper ion on azurin by analyzing both the apo and holo forms using temperature replica exchange molecular dynamics (REMD) data. We find a 14 K decrease in denaturation temperature for apo (406 K) azurin relative to that of holo (420 K), indicating a copper ion-induced additional thermal stability for holo azurin. The unfolding of apo azurin begins with the melting of α-helix and β-sheet V, similar to that of holo form. β-Sheets IV, VII, and VIII are comparatively more stable than other β-strands and melt at higher temperatures. Similar to holo azurin, the strong hydrophobic interactions among the apolar residues in the protein core is the key factor that renders high stability to apo protein as well. We construct free energy surfaces at different temperatures to capture the major conformations along the unfolding basins of the protein. Using contact maps from different basins we show the changes in the interaction between different residues along the unfolding pathway. Furthermore, we compare the Cα root-mean-square fluctuations (Cα-RMSF) and B-factor of all residues of apo and holo forms to understand the flexibility of different regions. The concerted displacement of α-helix and β-sheets V and VI from the protein core is another distinction we observe for apo compared to the holo form, where β-sheet VI was relatively stable.

Structure-based thermodynamics of ion selectivity (Mg2+versus Ca2+ and K+versus Na+) in the active site of the eukaryotic lariat group II intron from algae Pylaiella littoralis

Group II introns are metalloenzymes that can catalyze self-splicing. Recently, the crystal structures of the eukaryotic group IIB lariat intron from the brown algae Pylaiella littoralis have been reported for two intermediate states (pre-hydrolytic (2s) and post-hydrolytic) along the self-splicing pathway. Three characteristic metal-ion binding sites (M1 and M2 sites for catalytic Mg2+ ions, and K1 site for K+) in the catalytic pocket of the lariat intron have been identified and proposed to be crucial for self-splicing. Using the X-ray structures as a template, we quantitatively estimated the energetics of divalent (Mg2+versus Ca2+) and monovalent (K+versus Na+) ion selectivity and established a direct link between the energetics and structures of this lariat intron (bound to cognate and near-cognate metal ions). Molecular dynamics (MD) free energy simulations showed that the lariat intron was strongly selective between divalent metal ions. The strength of divalent metal-ion selectivity was noticeably high in the post-hydrolytic state (ΔΔG ≈ 20 kcal mol-1) relative to its pre-hydrolytic (2s) state (ΔΔG ≈ 13 kcal mol-1). Quantum chemical calculations ensured that the sign of the estimated divalent metal-ion selectivity was correct. The M1-binding pocket was less solvent-exposed in the case of the post-hydrolytic state relative to the pre-hydrolytic (2s) state, which boosted the metal-ion selectivity of the former. Surprisingly, in contrast to the bacterial linear group II intron, the lariat intron was found to be non-selective between monovalent ions (K+versus Na+). The interaction network in the first coordination shell of Ca2+ in the M1-binding pocket was different relative to Mg2+. Mg2+ substitution by Ca2+ resulted in the substitution of a single M1-RNA interaction by the M1-water interaction. In the pre-hydrolytic (2s) state, Ca2+ substitution completely disrupted the M1⋯5'-exon interaction; thus, the nature of the divalent metal ion is critical for catalysis. The interaction network in the M2 site was independent of the nature of the divalent metal ions (Mg2+ or Ca2+). The monovalent ion was loosely bound in the wet binding pocket (K1 site) of the lariat intron; thus, the substitution of K+ by Na+ could not significantly alter the free energy of the complex. The metal ion selectivity was dependent on the solvent accessibility of the metal-ion-binding-pocket, dry pocket enhanced the selectivity.

Modeling and dynamical analysis of the full-length structure of factor XII with zinc

Zinc (II), the second most abundant transition metal in blood, binds to the initiator of the contact pathway, factor XII (FXII). This binding induces conformational changes in the structure of FXII eventually leading to its activation. Despite many in vitro and in vivo studies on zinc-mediated activation of FXII, its molecular mechanism remains elusive mainly due to absence of a full-length structural model of FXII. To this end, this study investigated the role of zinc in the structure and dynamics of the full-length structure FXII that was obtained through molecular modeling. We have used four structural templates covering more than 70% of the FXII sequence and the remaining interconnecting regions were built by loop modeling. The resulting full-length structure of FXII contained disordered regions, but in comparison to the AlphaFold (AF) prediction, our full-length model represented a more realistic structure because of the disordered regions which were modeled to yield a more compact full-length structure in our model than the AF structure. Other than the disordered regions, our model and AF prediction were highly similar. The resulting full-length FXII structure was used to generate different systems representing the zinc-bound form (holo). Further to assess the contribution of the disulfide bridges, we also analyzed the apo and holo FXII structures with oxidized or reduced cysteine side-chains. Simulations suggested zinc binding conferred rigidity to the structure, particularly to the light chain of FXII. Zinc binding alone was sufficient to limit the backbone flexibility while 15 disulfide bonds, which were scattered throughout the structure, made a less significant contribution to the backbone rigidity. Altogether our results provide insights into the first realistic full-length structure of FXII focusing on the impact of structural zinc and disulfide bridges in the dynamics of this structure.

Review on the QM/MM Methodologies and Their Application to Metalloproteins

The multiscaling quantum mechanics/molecular mechanics (QM/MM) approach was introduced in 1976, while the extensive acceptance of this methodology started in the 1990s. The combination of QM/MM approach with molecular dynamics (MD) simulation, otherwise known as the QM/MM/MD approach, is a powerful and promising tool for the investigation of chemical reactions’ mechanism of complex molecular systems, drug delivery, properties of molecular devices, organic electronics, etc. In the present review, the main methodologies in the multiscaling approaches, i.e., density functional theory (DFT), semiempirical methodologies (SE), MD simulations, MM, and their new advances are discussed in short. Then, a review on calculations and reactions on metalloproteins is presented, where particular attention is given to nitrogenase that catalyzes the conversion of atmospheric nitrogen molecules N₂ into NH₃ through the process known as nitrogen fixation and the FeMo-cofactor.

Identification of allosteric fingerprints of alpha-synuclein aggregates in matrix metalloprotease-1 and substrate-specific virtual screening with single molecule insights

Alpha-synuclein (aSyn) has implications in pathological protein aggregations in neurodegeneration. Matrix metalloproteases (MMPs) are broad-spectrum proteases and cleave aSyn, leading to aggregation. Previous reports showed that allosteric communications between the two domains of MMP1 on collagen fibril and fibrin depend on substrates, activity, and ligands. This paper reports quantification of allostery using single molecule measurements of MMP1 dynamics on aSyn-induced aggregates by calculating Forster Resonance Energy Transfer (FRET) between two dyes attached to the catalytic and hemopexin domains of MMP1. The two domains of MMP1 prefer open conformations that are inhibited by a single point mutation E219Q of MMP1 and tetracycline, an MMP inhibitor. A two-state Poisson process describes the interdomain dynamics, where the two states and kinetic rates of interconversion between them are obtained from histograms and autocorrelations of FRET values. Since a crystal structure of aSyn-bound MMP1 is unavailable, binding poses were predicted by molecular docking of MMP1 with aSyn using ClusPro. MMP1 dynamics were simulated using predicted binding poses and compared with the experimental interdomain dynamics to identify an appropriate pose. The selected aSyn-MMP1 binding pose near aSyn residue K45 was simulated and analyzed to define conformational changes at the catalytic site. Allosteric residues in aSyn-bound MMP1 exhibiting strong correlations with the catalytic motif residues were compared with allosteric residues in free MMP1, and aSyn-specific residues were identified. The allosteric residues in aSyn-bound MMP1 are K281, T283, G292, G327, L328, E329, R337, F343, G345, N346, Y348, G353, Q354, D363, Y365, S366, S367, F368, P371, R372, V374, K375, A379, F391, A394, R399, M414, F419, V426, and C466. Shannon entropy was defined to quantify MMP1 dynamics. Virtual screening was performed against a site on selected aSyn-MMP1 binding poses, which showed that lead molecules differ between free MMP1 and substrate-bound MMP1. Also, identifying aSyn-specific allosteric residues in MMP1 enabled further selection of lead molecules. In other words, virtual screening needs to take substrates into account for potential substrate-specific control of MMP1 activity in the future. Molecular understanding of interactions between MMP1 and aSyn-induced aggregates may open up the possibility of degrading aggregates by targeting MMPs.

Molecular Insight into the High Thermal Stability of Metalloprotein Azurin

We investigate the events characterizing the steps of the unfolding pathway of blue copper metalloprotein azurin using replica exchange molecular dynamics (REMD). Our studies show that the unfolding of azurin begins with the melting of α-helix and β-sheets II and V. This is followed by the melting of other β-sheets and the exposure of hydrophobic protein core to the solvent, resulting in disruptions of its tertiary structure. Free energy surfaces constructed at different temperatures portray different basins that signify the stability of different melted structures in the unfolding process. The contact maps at different temperatures reveal that the strong hydrophobic interaction within the core of the protein is the vital force that renders high stability to this protein. Analysis of the individual β-sheets by looking into their amino acid sequence shows that β-sheets with charged side chains on the surface melt fast compared to others. The β-barrel of azurin is able to dynamically rearrange, and it helps the protein to preserve its hydrophobic core, holding back the native topology from melting fast. B-factor analysis shows that residues of β-sheets III, IV, and VII deviate less from their initial structure at the transition temperature.

Divalent-Metal-Ion Selectivity of the CRISPR-Cas System-Associated Cas1 Protein: Insights from Classical Molecular Dynamics Simulations and Electronic Structure Calculations

CRISPR-associated protein 1 (Cas1) is a universally conserved essential metalloenzyme of the clustered regularly interspaced short palindromic repeat (CRISPR) immune system of prokaryotes (bacteria, archaea) that can cut and integrate a part of viral DNA to its host genome with the help of other proteins. The integrated DNA acts as a memory of viral infection, which can be transcribed to RNA and stop future infection by recognition (based on the RNA/DNA complementarity principle) followed by protein-mediated degradation of the viral DNA. It has been proposed that the presence of a single manganese (Mn²⁺) ion in a conserved divalent-metal-ion binding pocket (key residues: E190, H254, D265, D268) of Cas1 is crucial for its function. Cas1-mediated DNA degradation was proposed to be hindered by metal substitution, metal chelation, or mutation of the binding pocket residues. Cas1 is active toward dsDNA degradation with both Mn²⁺ and Mg²⁺. X-ray structures of Cas1 revealed an intricate atomic interaction network of the divalent-metal-ion binding pocket and opened up the possibility of modeling related metal ions (viz., Mg²⁺, Ca²⁺) in the binding pocket of wild-type (WT) and mutated Cas1 proteins for computational analysis, which includes (1) quantitative estimation of the energetics of the divalent-metal-ion preference and (2) exploring the structural and dynamical aspects of the protein in response to divalent-metal-ion substitution or amino acid mutation. Using the X-ray structure of the Cas1 protein from Pseudomonas aeruginosa as a template (PDB 3GOD), we performed (∼2.23 μs) classical molecular dynamics (MD) simulations to compare structural and dynamical differences between Mg²⁺- and Ca²⁺-bound binding pockets of wild-type (WT) and mutant (E190A, H254A, D265A, D268A) Cas1. Furthermore, reduced binding pocket models were generated from X-ray and molecular dynamics (MD) trajectories, and the resulting structures were subjected to quantum chemical calculations. Results suggest that Cas1 prefers Mg²⁺ binding relative to Ca²⁺ and the preference is the strongest for WT and the weakest for the D268A mutant. Quantum chemical calculations indicate that Mn²⁺ is the most preferred relative to both Mg²⁺ and Ca²⁺ in the wild-type and mutant Cas1. Substitution of Mg²⁺ by Ca²⁺ does not alter the interaction network between Cas1 and the divalent metal ion but increases the wetness of the binding pocket by introducing a single water molecule in the first coordination shell of the latter. The strength of metal-ion preference (Mg²⁺ versus Ca²⁺) seems to be dependent on the solvent accessibility of the divalent-metal-ion binding pocket, strongest for wild-type Cas1 (in which the metal-ion binding pocket is dry, which includes two water molecules) and the weakest for the D268A mutant (in which the metal-ion binding pocket is wet, which includes four water molecules).

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.

Finding High-Quality Metal Ion-Centric Regions Across the Worldwide Protein Data Bank

As the number of macromolecular structures in the worldwide Protein Data Bank (wwPDB) continues to grow rapidly, more attention is being paid to the quality of its data, especially for use in aggregated structural and dynamics analyses. In this study, we systematically analyzed 3.5 Å regions around all metal ions across all PDB entries with supporting electron density maps available from the PDB in Europe. All resulting metal ion-centric regions were evaluated with respect to four quality-control criteria involving electron density resolution, atom occupancy, symmetry atom exclusion, and regional electron density discrepancy. The resulting list of metal binding sites passing all four criteria possess high regional structural quality and should be beneficial to a wide variety of downstream analyses. This study demonstrates an approach for the pan-PDB evaluation of metal binding site structural quality with respect to underlying X-ray crystallographic experimental data represented in the available electron density maps of proteins. For non-crystallographers in particular, we hope to change the focus and discussion of structural quality from a global evaluation to a regional evaluation, since all structural entries in the wwPDB appear to have both regions of high and low structural quality.

Computational approaches for deciphering the equilibrium and kinetic properties of iron transport proteins

With the advances in three-dimensional structure determination techniques, high quality structures of the iron transport proteins transferrin and the bacterial ferric binding protein (FbpA) have been deposited in the past decade. These are proteins of relatively large size, and developments in hardware and software have only recently made it possible to study their dynamics using standard computational resources. We review computational techniques towards understanding the equilibrium and kinetic properties of iron transport proteins under different environmental conditions. At the level of detail that requires quantum chemical treatments, the octahedral geometry around iron has been scrutinized and it has been established that the iron coordinating tyrosines are in an unusual deprotonated state. At the atomistic level, both the N-lobe and the full bilobal structure of transferrin have been studied under varying conditions of pH, ionic strength and binding of other metal ions by molecular dynamics (MD) simulations. These studies have allowed questions to be answered, among others, on the function of second shell residues in iron release, the role of synergistic anions in preparing the active site for iron binding, and the differences between the kinetics of the N- and the C-lobe. MD simulations on FbpA have led to the detailed observation of the binding kinetics of phosphate to the apo form, and to the conformational preferences of the holo form under conditions mimicking the environmental niches provided by the periplasmic space. To study the dynamics of these proteins with their receptors, one must resort to coarse-grained methodologies, since these systems are prohibitively large for atomistic simulations. A study of the complex of human transferrin (hTf) with its pathogenic receptor by such methods has revealed a potential mechanistic explanation for the defense mechanism that arises in evolutionary warfare. Meanwhile, the motions in the transferrin receptor bound hTf have been shown to disfavor apo hTf dissociation, explaining why the two proteins remain in complex during the recycling process from the endosome to the cell surface. Open problems and possible technological applications related to metal ion binding-release in iron transport proteins that may be handled by hybrid use of quantum mechanical, MD and coarse-grained approaches are discussed.

Challenges and advances in the computational modeling of biological phosphate hydrolysis

Phosphate ester hydrolysis is fundamental to many life processes, and has been the topic of substantial experimental and computational research effort. However, even the simplest of phosphate esters can be hydrolyzed through multiple possible pathways that can be difficult to distinguish between, either experimentally, or computationally. Therefore, the mechanisms of both the enzymatic and non-enzymatic reactions have been historically controversial. In the present contribution, we highlight a number of technical issues involved in reliably modeling these computationally challenging reactions, as well as proposing potential solutions. We also showcase examples of our own work in this area, discussing both the non-enzymatic reaction in aqueous solution, as well as insights obtained from the computational modeling of organophosphate hydrolysis and catalytic promiscuity amongst enzymes that catalyze phosphoryl transfer.

Selective cysteine modification of metal-free human metallothionein 1a and its isolated domain fragments: Solution structural properties revealed via ESI-MS

Human metallothionein 1a, a protein with two cysteine-rich metal-binding domains (α with 11 Cys and β with 9), was analyzed in its metal-free form by selective, covalent Cys modification coupled with ESI-MS. The modification profiles of the isolated β- and α-fragments reacted with p-benzoquinone (Bq), N-ethylmalemide (NEM) and iodoacetamide (IAM) were compared with the full length protein using ESI-mass spectral data to follow the reaction pathway. Under denaturing conditions at low pH, the reaction profile with each modifier followed pathways that resulted in stochastic, Normal distributions of species whose maxima was equal to the mol. eq. of modifier added. Our interpretation of modification at this pH is that reaction with the cysteines is unimpeded when the full protein or those of its isolated domains are denatured. At neutral pH, where the protein is expected to be folded in a more compact structure, there is a difference in the larger Bq and NEM modification, whose reaction profiles indicate a cooperative pattern. The reaction profile with IAM under native conditions follows a similar stochastic distribution as at low pH, suggesting that this modifier is small enough to access the cysteines unimpeded by the compact structure. The data emphasize the utility of residue modification coupled with electrospray ionization mass spectrometry for the study of protein structure.

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.

Developing multisite empirical force field models for Pt(II) and cisplatin

We have developed empirical force field parameters for Pt(II) and cisplatin. Two force field frameworks were used-modified OPLS-AA and our second-order polarizable POSSIM. A seven-site model was used for the Pt(II) ion. The goal was to create transferable parameter sets compatible with the force field models for proteins and general organic compounds. A number of properties of the Pt(II) ion and its coordination compounds have been considered, including geometries and energies of the complexes, hydration free energy, and radial distribution functions in water. Comparison has been made with experimental and quantum mechanical results. We have demonstrated that both versions are generally capable of reproducing key properties of the system, but the second-order polarizable POSSIM formalism permits more accurate quantitative results to be obtained. For example, the energy of formation of cisplatin as calculated with the modified OPLS-AA exhibited an error of 9.9%, while the POSSIM error for the same quantity was 6.2%. The produced parameter sets are transferable and suitable to be used in protein-metal binding simulations in which position or even coordination of the ion does not have to be constrained using preexisting knowledge. © 2016 Wiley Periodicals, Inc.

Zn Coordination Chemistry: Development of Benchmark Suites for Geometries, Dipole Moments, and Bond Dissociation Energies and Their Use To Test and Validate Density Functionals and Molecular Orbital Theory

We present nonrelativistic and relativistic benchmark databases (obtained by coupled cluster calculations) of 10 Zn-ligand bond distances, 8 dipole moments, and 12 bond dissociation energies in Zn coordination compounds with O, S, NH3, H2O, OH, SCH3, and H ligands. These are used to test the predictions of 39 density functionals, Hartree-Fock theory, and seven more approximate molecular orbital theories. In the nonrelativisitic case, the M05-2X, B97-2, and mPW1PW functionals emerge as the most accurate ones for this test data, with unitless balanced mean unsigned errors (BMUEs) of 0.33, 0.38, and 0.43, respectively. The best local functionals (i.e., functionals with no Hartree-Fock exchange) are M06-L and τ-HCTH with BMUEs of 0.54 and 0.60, respectively. The popular B3LYP functional has a BMUE of 0.51, only slightly better than the value of 0.54 for the best local functional, which is less expensive. Hartree-Fock theory itself has a BMUE of 1.22. The M05-2X functional has a mean unsigned error of 0.008 Å for bond lengths, 0.19 D for dipole moments, and 4.30 kcal/mol for bond energies. The X3LYP functional has a smaller mean unsigned error (0.007 Å) for bond lengths but has mean unsigned errors of 0.43 D for dipole moments and 5.6 kcal/mol for bond energies. The M06-2X functional has a smaller mean unsigned error (3.3 kcal/mol) for bond energies but has mean unsigned errors of 0.017 Å for bond lengths and 0.37 D for dipole moments. The best of the semiempirical molecular orbital theories are PM3 and PM6, with BMUEs of 1.96 and 2.02, respectively. The ten most accurate functionals from the nonrelativistic benchmark analysis are then tested in relativistic calculations against new benchmarks obtained with coupled-cluster calculations and a relativistic effective core potential, resulting in M05-2X (BMUE = 0.895), PW6B95 (BMUE = 0.90), and B97-2 (BMUE = 0.93) as the top three functionals. We find significant relativistic effects (∼0.01 Å in bond lengths, ∼0.2 D in dipole moments, and ∼4 kcal/mol in Zn-ligand bond energies) that cannot be neglected for accurate modeling, but the same density functionals that do well in all-electron nonrelativistic calculations do well with relativistic effective core potentials. Although most tests are carried out with augmented polarized triple-ζ basis sets, we also carried out some tests with an augmented polarized double-ζ basis set, and we found, on average, that with the smaller basis set DFT has no loss in accuracy for dipole moments and only ∼10% less accurate bond lengths.

Insights into the key interactions between human protein phosphatase 5 and cantharidin using molecular dynamics and site-directed mutagenesis bioassays

Serine/threonine protein phosphatase 5 (PP5) is a promising novel target for anticancer therapies. This work aims to uncover the key interactions at the atomic level between PP5 and three inhibitors (cantharidin, norcantharidin and endothall). We found that, unlike previous report, Arg 100 contributes less to PP5-inhibitor binding, and the residues His 69, Asn 128, His 129, Arg 225, His 252 and Arg 250 are of importance to PP5-inhibitor binding. The hydrophobic interactions established between the residues Val 254, Phe 271 and Tyr 276, especially Glu 253, are very important to enhance the inhibitive interaction. We suggested that, to increase the inhibitory activity, the interactions of inhibitor with three negatively charged unfavorable interaction residues, Asp 99, Glu 130 and Asp 213, should be avoided. However, the interactions of inhibitor with favorable interaction residue Arg 250 could enhance the inhibitory activity. The Manganese ion 2 (MN2) unfavorably contribute to the total interaction free energies. The coordination between MN2 and chemical group of inhibitor should be eliminated. This work provides insight into how cantharidin and its analogs bind to PP5c at the atomic level and will facilitate modification of cantharidin-like chemicals to rationally develop more specific and less cytotoxic anti-cancer drugs.

Integrating atomistic molecular dynamics simulations, experiments, and network analysis to study protein dynamics: strength in unity

In the last years, we have been observing remarkable improvements in the field of protein dynamics. Indeed, we can now study protein dynamics in atomistic details over several timescales with a rich portfolio of experimental and computational techniques. On one side, this provides us with the possibility to validate simulation methods and physical models against a broad range of experimental observables. On the other side, it also allows a complementary and comprehensive view on protein structure and dynamics. What is needed now is a better understanding of the link between the dynamic properties that we observe and the functional properties of these important cellular machines. To make progresses in this direction, we need to improve the physical models used to describe proteins and solvent in molecular dynamics, as well as to strengthen the integration of experiments and simulations to overcome their own limitations. Moreover, now that we have the means to study protein dynamics in great details, we need new tools to understand the information embedded in the protein ensembles and in their dynamic signature. With this aim in mind, we should enrich the current tools for analysis of biomolecular simulations with attention to the effects that can be propagated over long distances and are often associated to important biological functions. In this context, approaches inspired by network analysis can make an important contribution to the analysis of molecular dynamics simulations.

Importance of polarization effect in the study of metalloproteins: application of polarized protein specific charge scheme in predicting the reduction potential of azurin

Molecular dynamics (MD) simulation is commonly used in the study of protein dynamics, and in recent years, the extension of MD simulation to the study of metalloproteins is gaining much interest. Choice of force field is crucial in MD studies, and the inclusion of metal centers complicates the process of accurately describing the electrostatic environment that surrounds the redox centre. Herein, we would like to explore the importance of including electrostatic contribution from both protein and solvent in the study of metalloproteins. MD simulations with the implementation of thermodynamic integration will be conducted to model the reduction process of azurin from Pseudomonas aeruginosa. Three charge schemes will be used to derive the partial charges of azurin. These charge schemes differ in terms of the amount of immediate environment, respective to copper, considered during charge fitting, which ranges from the inclusion of copper and residues in the first coordination sphere during density functional theory charge fitting to the comprehensive inclusion of protein and solvent effect surrounding the metal centre using polarized protein-specific charge scheme. From the simulations conducted, the relative reduction potential of the mutated azurins respective to that of wild-type azurin (ΔEcal) were calculated and compared with experimental values. The ΔEcal approached experimental value with increasing consideration of environmental effect hence substantiating the importance of polarization effect in the study of metalloproteins. This study also attests the practicality of polarized protein-specific charge as a computational tool capable of incorporating both protein environment and solvent effect into MD simulations.

Modeling catalytic promiscuity in the alkaline phosphatase superfamily

In recent years, it has become increasingly clear that promiscuity plays a key role in the evolution of new enzyme function. This finding has helped to elucidate fundamental aspects of molecular evolution. While there has been extensive experimental work on enzyme promiscuity, computational modeling of the chemical details of such promiscuity has traditionally fallen behind the advances in experimental studies, not least due to the nearly prohibitive computational cost involved in examining multiple substrates with multiple potential mechanisms and binding modes in atomic detail with a reasonable degree of accuracy. However, recent advances in both computational methodologies and power have allowed us to reach a stage in the field where we can start to overcome this problem, and molecular simulations can now provide accurate and efficient descriptions of complex biological systems with substantially less computational cost. This has led to significant advances in our understanding of enzyme function and evolution in a broader sense. Here, we will discuss currently available computational approaches that can allow us to probe the underlying molecular basis for enzyme specificity and selectivity, discussing the inherent strengths and weaknesses of each approach. As a case study, we will discuss recent computational work on different members of the alkaline phosphatase superfamily (AP) using a range of different approaches, showing the complementary insights they have provided. We have selected this particular superfamily, as it poses a number of significant challenges for theory, ranging from the complexity of the actual reaction mechanisms involved to the reliable modeling of the catalytic metal centers, as well as the very large system sizes. We will demonstrate that, through current advances in methodologies, computational tools can provide significant insight into the molecular basis for catalytic promiscuity, and, therefore, in turn, the mechanisms of protein functional evolution.

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.

Carbonic anhydrase binding site parameterization in OPLS-AA force field

The parameterization of carbonic anhydrase binding site in OPLS-AA force field was performed using quantum chemistry calculations. Both OH2 and OH(-) forms of the binding site were considered, showing important differences in terms of atomic partial charges. Three different parameterization protocols were used, and the results obtained highlighted the importance of including an extended binding site in the charge calculation. The force field parameters were subsequently validated using standard molecular dynamics simulations. The results presented in this work should greatly facilitate the use of molecular dynamics simulations for studying the carbonic anhydrase, and more generally, the metallo-enzymes.

The reversible opening of water channels in cytochrome c modulates the heme iron reduction potential

Dynamic protein-solvent interactions are fundamental for life processes, but their investigation is still experimentally very demanding. Molecular dynamics simulations up to hundreds of nanoseconds can bring to light unexpected events even for extensively studied biomolecules. This paper reports a combined computational/experimental approach that reveals the reversible opening of two distinct fluctuating cavities in Saccharomyces cerevisiae iso-1-cytochrome c. Both channels allow water access to the heme center. By means of a mixed quantum mechanics/molecular dynamics (QM/MD) theoretical approach, the perturbed matrix method (PMM), that allows to reach long simulation times, changes in the reduction potential of the heme Fe(3+)/Fe(2+) couple induced by the opening of each cavity are calculated. Shifts of the reduction potential upon changes in the hydration of the heme propionates are observed. These variations are relatively small but significant and could therefore represent a tool developed by cytochrome c for the solvent driven, fine-tuning of its redox functionality.

Structures and free energy landscapes of aqueous zinc(II)-bound amyloid-β(1-40) and zinc(II)-bound amyloid-β(1-42) with dynamics

Binding of divalent metal ions with intrinsically disordered fibrillogenic proteins, such as amyloid-β (Aβ), influences the aggregation process and the severity of neurodegenerative diseases. The Aβ monomers and oligomers are the building blocks of the aggregates. In this work, we report the structures and free energy landscapes of the monomeric zinc(II)-bound Aβ40 (Zn:Aβ40) and zinc(II)-bound Aβ42 (Zn:Aβ42) intrinsically disordered fibrillogenic metallopeptides in an aqueous solution by utilizing an approach that employs first principles calculations and parallel tempering molecular dynamics simulations. The structural and thermodynamic properties, including the secondary and tertiary structures and conformational Gibbs free energies of these intrinsically disordered metallopeptide alloforms, are presented. The results show distinct differing characteristics for these metallopeptides. For example, prominent β-sheet formation in the N-terminal region (Asp1, Arg5, and Tyr10) of Zn:Aβ40 is significantly decreased or lacking in Zn:Aβ42. Our findings indicate that blocking multiple reactive residues forming abundant β-sheet structure located in the central hydrophobic core and C-terminal regions of Zn:Aβ42 via antibodies or small organic molecules might help to reduce the aggregation of Zn(II)-bound Aβ42. Furthermore, we find that helix formation increases but β-sheet formation decreases in the C-terminal region upon Zn(II) binding to Aβ. This depressed β-sheet formation in the C-terminal region (Gly33-Gly38) in monomeric Zn:Aβ42 might be linked to the formation of amorphous instead of fibrillar aggregates of Zn:Aβ42.

Importance of electrostatic polarizability in calculating cysteine acidity constants and copper(I) binding energy of Bacillus subtilis CopZ

CopZ is a copper chaperone from Bacillus subtilis. It is an important part of Cu(I) trafficking. We have calculated pK(a) values for the CXXC motif of this protein, which is responsible for the Cu(I) binding, and the Cu(I) binding constants. Polarizable and fixed-charges formalisms were used, and solvation parameters for the both models have been refitted. We had to partially redevelop parameters for the protonated and deprotonated cysteine residues. We have discovered that the polarizable force field (PFF) is qualitatively superior and allows a uniformly better level of energetic results. The PFF pK(a) values for cysteine are within about 0.8-2.8 pH units of the experimental data, while the fixed-charges OPLS formalism yields errors of up to tens of units. The PFF magnitude of the copper binding energy is about 10 kcal/mol or 50% higher than the experimental value, while the using the refitted OPLS parameters leads to an overall positive binding energy, thus predicting no thermodynamically stable complex. At the same time, the agreement of the polarizable S···Cu(I) distances with the experimental results is within 0.08 Å range, and the nonpolarizable calculations lead to an error of about 0.4 Å. Moreover, the accuracy of the PFF has been achieved without any explicit fitting to either pK(a) or CopZ···Cu(I) binding energies. We believe that this makes our polarizable technique a choice method in reproducing protein-copper binding and further supports the notion that explicit treatment of electrostatic polarization is crucial in many biologically relevant studies, especially ion binding and transport.

Electrostatic polarization is crucial in reproducing Cu(I) interaction energies and hydration

We have explored the suitability of fixed-charges and polarizable force fields for modeling interactions of the monovalent Cu(I) ion. Parameters for this ion have been tested and refitted within the fixed-charges OPLS-AA and polarizable force field (PFF) frameworks. While this ion plays an important role in many protein interactions, the attention to it in developing empirical force fields is limited. Our PFF parameters for the copper ion worked very well for the Cu(I) interactions with water, while both the original OPLS2005 and our refitted OPLS versions moderately underestimated the copper-water interaction energy. However, the greatest problem in using the nonpolarizable fixed-charges OPLS force field was observed while calculating interaction energies and distances for Cu(I)-benzene complexes. The OPLS2005 model underestimates the interaction energy by a factor of 4. Refitting the OPLS parameters reduced this underestimation to a factor of 2.2-2.4, but only at a cost of distorting the complex geometry. At the same time, the polarizable calculations had an error of about 4%. Moreover, we then used the PFF and nonpolarizable refitted OPLS models for finding free energy of hydration for copper ion via molecular dynamics simulations. While the OPLS calculations lead to a 22% error in the solvation energy, the PFF result was off by only 1.8%. This was achieved with no refitting of the parameters but simply by employing the model developed for the Cu(I) interaction with a single water molecule. We believe that the presented results not only lead to a conclusion about a qualitatively greater suitability of polarizable force fields for simulating molecular interactions with ions but also attest to the excellent level of transferability of PFF parameters.

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.

Structure and dynamics of Cu(I) binding in copper chaperones Atox1 and CopZ: a computer simulation study

Copper chaperones deliver reduced copper (i.e., Cu(I)) to metal-binding domains of P-type ATPases in the cytoplasm of a range of organisms. Both chaperones and target domains have a ferredoxin-like fold and metal-binding motifs involving two Cys residues. Here, we investigated the Cu-binding geometry and structural dynamics of two homologous Cu(I) chaperones, Homo sapiens Atox1 and Bacillus subtilis CopZ, using a combination of quantum mechanical-molecular mechanics (QM-MM) and classical molecular dynamics (MD) methods. Our QM-MM optimized geometries for the holo- proteins suggested that Cu(I) in Atox1 favors a linear Cys(S)-Cu-Cys(S) arrangement but that this angle is close to 150 degrees in CopZ. Classical MD simulations suggest that both Atox1 and CopZ apo- forms have an increased conformational flexibility as compared to the respective holo- forms. This difference is most pronounced in CopZ and correlates with a lower in vitro thermal stability. Both average fluctuation (i.e., rmsd) and radius of gyration data demonstrate that the effects of Cu(I) coordination extend throughout the proteins. Distinct deviations between the two homologues were found in protein-solvent interactions, entropy of Cu(I) binding, and apo-protein Cys-Cys distance distributions. Our in silico results provide new insights into copper chaperone behavior with direct implications for copper transport mechanisms in vivo.

Structural prediction of the beta-domain of metallothionein-3 by molecular dynamics simulation

The beta-domain of metallothionein-3 (MT3) has been reported to be crucial to the neuron growth inhibitory bioactivity. Little detailed three-dimensional structural information is available to present a reliable basis for elucidation on structure-property-function relationships of this unique protein by experimental techniques. So, molecular dynamics simulation is adopted to study the structure of beta-domain of MT3. In this article, a 3D structural model of beta-domain of MT3 was generated. The molecular simulations provide detailed protein structural information of MT3. As compared with MT2, we found a characteristic conformation formed in the fragment (residue 1-13) at the N-terminus of MT3 owing to the constraint induced by 5TCPCP9, in which Pro7 and Pro9 residues are on the same side of the protein, both facing outward and the two 5-member rings of prolines are arranged almost in parallel, while Thr5 is on the opposite side. Thr5 in MT3 is also found to make the first four residues relatively far from the fragment (residue 23-26) as compared with MT2. The simulated structure of beta-domain of MT3 is looser than that of MT2. The higher energy of MT3 than that of MT2 calculated supports these conclusions. Simulation on the four isomer arising from the cis- or trans-configuration of 6CPCP9 show that the trans-/trans-isomer is energetic favorable. The partially unfolding structure of beta-domain of MT3 is also simulated and the results show the influence of 6CPCP9 sequence on the correct folding of this domain. The correlations between the bioactivity of MT3 and the simulated structure as well as the folding of beta-domain of MT3 are discussed based on our simulation and previous results.

Comparative assessment of theoretical methods for the determination of geometrical properties in biological zinc complexes

In the present study, we have compared the performance of the density functional theory (DFT) functionals B1B95, B3LYP, B97-2, BP86, and BPW91 with MP2 for geometry determination in biological mononuclear Zn complexes. A total of 15 different basis sets, of rather diverse complexity, were tested, several which included also three different types of common effective-core potentials: Los Alamos, Steven-Basch-Krauss, and Stuttgart-Dresden. In addition, the ability to describe mononuclear Zn biological systems using relatively simple models of the metal coordination sphere, comprising only the metal atom and a simplified representation of the ligands at the first coordination sphere, starting from a set of high-resolution X-ray crystallographic structures, is evaluated for 90 combinations of method/basis set. The results show that the use of such models allows for a relatively accurate description of the Zn-ligand bond lengths, although failing to correctly represent the topology of the metal coordination sphere (namely, the angles involving the metal atom) if constraints at the Calpha atoms are not considered. Globally, B3LYP had the best average performance in the test, closely followed by MP2, whereas B1B95 was the least accurate method. The study also points out B3LYP/CEP-121G and B3LYP/SDD, which use, respectively, the Steven-Basch-Krauss and the Stuttgart-Dresden effective-core potentials, as the best compromise between accuracy and CPU time for the geometrical characterization of metal-ligand bond lengths in Zn biological systems.

Computational protocol for predicting the binding affinities of zinc containing metalloprotein-ligand complexes

Zinc is one of the most important metal ions found in proteins performing specific functions associated with life processes. Coordination geometry of the zinc ion in the active site of the metalloprotein-ligand complexes poses a challenge in determining ligand binding affinities accurately in structure-based drug design. We report here an all atom force field based computational protocol for estimating rapidly the binding affinities of zinc containing metalloprotein-ligand complexes, considering electrostatics, van der Waals, hydrophobicity, and loss in conformational entropy of protein side chains upon ligand binding along with a nonbonded approach to model the interactions of the zinc ion with all the other atoms of the complex. We examined the sensitivity of the binding affinity predictions to the choice of Lennard-Jones parameters, partial atomic charges, and dielectric treatments adopted for system preparation and scoring. The highest correlation obtained was R2 = 0.77 (r = 0.88) for the predicted binding affinity against the experiment on a heterogenous dataset of 90 zinc containing metalloprotein-ligand complexes consisting of five unique protein targets. Model validation and parameter analysis studies underscore the robustness and predictive ability of the scoring function. The high correlation obtained suggests the potential applicability of the methodology in designing novel ligands for zinc-metalloproteins. The scoring function has been web enabled for free access at www.scfbio-iitd.res.in/software/drugdesign/bapplz.jsp as BAPPL-Z server (Binding Affinity Prediction of Protein-Ligand complexes containing Zinc metal ions).

The Characterization and Role of Zinc Binding in Yeast Cox4

Yeast Cox4 is a zinc binding subunit of cytochrome c oxidase. Cox4 is the only cofactor-containing subunit that is not directly part of the catalytic core of the enzyme located in the mitochondrial inner membrane. The Zn(II) site is shown to be distinct from the bovine ortholog, as it results from the x-ray structure of the entire cytochrome c oxidase in having a single histidyl residue and three conserved cysteines residues in the coordination sphere. Substitutions at the Cys ligand positions result in non-functional Cox4 proteins that fail to lead to cytochrome oxidase assembly. Limited function exists in His-119 mutants when overexpressed. Zn(II) binding in Cox4 is, therefore, important for the stability of the complex. The solution structure of yeast Cox4 elucidated by multidimensional NMR reveals a C-terminal globular domain consisting of two beta sheets analogous to the bovine ortholog except the loop containing the coordinating His in the yeast protein and the fourth Cys in the bovine protein are in different positions in the two structures. The conformation of this loop is dictated by the different sequence position of the fourth coordinating zinc ligand. The Zn(II) ion is buried within the domain, consistent with its role in structural stability. Potential functions of this matrix-facing subunit are discussed.

Modeling the charge distribution at metal sites in proteins for molecular dynamics simulations

Almost half of the proteome of living organisms is constituted of metalloproteins. Unfortunately, the ability of the current generation of molecular dynamics pairwise-additive forcefields to properly describe metal pockets is severely lacking due to the intrinsic difficulty of handling polarization and charge transfer contributions. In order to improve the description of metalloproteins, a simple reparameterization strategy is proposed herein that does not involve artificial constraints. Specifically, a non-bonded quantum mechanical-based model is used to capture the mean polarization and charge transfer contributions to the interatomic forces within the metal site. The present approach is demonstrated to provide enough accuracy to maintain the integrity of the metal pocket for a variety of metalloproteins during extended (multi-nanosecond) molecular dynamics simulations. The method enables the sampling of small conformational changes and the relaxation of local frustrations in NMR structures.