Int J App Pharm, Vol 16, Issue 5, 2024, 90-98Original Article

METALLOPROTEIN PARAMETERS IN MOLECULAR DYNAMICS SIMULATION FOR AMBER, CHARMM, GROMACS, AND NAMD: A SYSTEMATIC REVIEW

PURNAWAN PONTANA PUTRA1*, NAJMIATUL FITRIA2, AIYI ASNAWI3, AKMAL DJAMAAN1

1Department of Pharmaceutical Chemistry, Faculty of Pharmacy, Universitas Andalas, Padang-25163, Indonesia. 2Department of Pharmacology and Clinical Pharmacy, Faculty of Pharmacy, Universitas Andalas, Padang-25163, Indonesia. 3Faculty of Pharmacy, Bhakti Kencana University, Bandung-40614, Indonesia
*Corresponding author: Purnawan Pontana Putra; *Email: purnawanpp@phar.unand.ac.id

Received: 20 May 2024, Revised and Accepted: 12 Jul 2024

ABSTRACT

Objective: The selection of appropriate metal parameters for molecular dynamics simulations is a significant challenge. Therefore, this review aims to provide in-depth insights valuable for the optimization of parameter selection in the context of chemical simulations.

Methods: A total of 550 scientific articles were collected from pubmed and science direct databases from 2009 to 2024, resulting in the inclusion of 60 full studies for review. The selection process of preferred reporting items for systematic reviews and meta-analyses (PRISMA) was utilized, enabling the conduction of an initial screening of articles by use of the Rayyan web-based application.

Results: This study found that the modeling and parameterization of metal proteins were categorized into bonded and non-bonded models. The Bonded Model incorporates MCPB, a Python-based software that facilitates parameter construction for over 80 metal ions and force fields in molecular dynamics simulations. The non-bonded model evaluates metals in proteins, such as zinc, nickel, magnesium, cobalt, iron, and cadmium by using AMBER force field and the Seminario method. The 12-6 lennard-Jones (LJ) non-bonded model is suitable for divalent, trivalent, and tetravalent metals, with Zinc parameters being compared for accuracy. Additionally, the force fields suitable for modeling unbound metal proteins include AMBER FF19SB, FF14SB, ff9X, CHARMM36, CHARMM22, CHARMM27, and CHARMM-Metal.

Conclusion: This study found that the modeling and parameterization of metal proteins were categorized into bonded and non-bonded models. molecular dynamics (MD) simulations can be conducted using various methods, such as classical molecular dynamics, umbrella sampling, quantum mechanics-discrete molecular dynamics (QM/DMD), stochastic boundary molecular dynamics (SBMD), steered molecular dynamics (SMD), gaussian accelerated molecular dynamics (GaMD) and random acceleration molecular dynamics (RAMD).

Keywords: Metalloprotein, Parameters, Molecular dynamics, Systematic review

© 2024 The Authors. Published by Innovare Academic Sciences Pvt Ltd. This is an open access article under the CC BY license (https://creativecommons.org/licenses/by/4.0/)
DOI: https://dx.doi.org/10.22159/ijap.2024v16i5.51513 Journal homepage: https://innovareacademics.in/journals/index.php/ijap

INTRODUCTION

The role of metals in various reactions is essential, thereby attracting significant attention to modeling atomic interactions. In understanding how metals interact with other molecules, molecular dynamics simulations are key. However, using the right parameters is crucial to accurately show how atoms move. In a Bonded Model, metals in proteins are typically placed and linked to ligands, whereas in a non-bonded model, metals are free from ligands [1]. Previous studies developed computational methods to find genes encoding metalloproteins and to identify metal-binding sites or patterns in proteins, along with introducing related databases. Moreover, recent research discusses advancements in comparing how important metals are processed differently in prokaryotes and eukaryotes, revealing diverse evolutionary patterns in various metalloprotein families and metalloproteomes [2].

Several challenges in molecular dynamics simulations of metal proteins can be analyzed through a de novo method [3]. Additionally, parameterization for metals is necessary before conducting molecular dynamics simulations using software such as AMBER, GROMACS, CHARMM, and NAMD. The most frequently used parameterization method is Quantum Quantum Mechanics/Molecular Mechanics (QM/MM). This approach preserves the realistic description of the protein environment assessed at more computationally tractable molecular mechanics level. Previous studies focused on modeling metal ions with Quantum mechanical methods and force field models [4]. Parameterizing metals is crucial as it enables the calculation of protein-ligand or organometallic compound affinities, which is essential for understanding their interactions in biological and chemical contexts [5, 6]. Moreover, the merging of quantum mechanics (QM) and molecular mechanics (MM) in simulation offers an ideal means of providing a true quantum-chemical description of complicated structures comprising metal coordination spheres.

Different types of metals require different analysis and accuracy from the available collection of parameters. Greater accuracy usually implies more complex parameters and longer consideration times [7]. Various metal parameters may also be more suitable for certain simulation purposes, such as the study of the structure of metal-protein complexes, electron transport, or certain chemical reactions [8]. The molecular dynamics simulation parameters are optimized for several evaluation metrics using a collection of computational methods in pursuit of the best balance between accuracy and efficiency. More importantly, proper choice of the parameters-quantum and classical modeling strategy is necessary to generate representative and relevant results [9].

Standardized guidelines are increasingly needed for the selection and adjustment of metal parameters in molecular dynamics simulations performed with AMBER, CHARMM, GROMACS, and NAMD software, necessitating the exploration of previous studies. Therefore, this review aims to bridge the knowledge gap regarding the application of metal parameters in molecular dynamics simulations using the four software programs mentioned. Valuable comprehensive insights are expected to be provided for the optimization of metal parameter selection in the context of chemical simulations.

MATERIALS AND METHODS

Search strategy

Approximately 413 and 137 articles (totaling 550) were respectively collected from Pubmed and ScienceDirect databases, covering the period of 2009 to 2024. The collection process incorporated inclusion and exclusion criteria, filtering, data extraction, and assessment methods [10]. Subsequently, screening was conducted using the preferred reporting items for systematic reviews and meta-analyses (PRISMA) method [11], and articles were selected with the web-based application known as Rayyan [12]. Literature identification was accomplished through truncated searches using combinations of (("molecular dynamics simulation"[MeSH Terms] OR molecular dynamic [Text Word]) AND "Molecular Dynamics Simulation"[Mesh]) AND (("metalloproteins"[MeSH Terms] OR metalloprotein[Text Word]) AND "Metalloproteins"[Mesh]).

Fig. 1: PRISMA flowchart depicting the screening and selection process of the published literature

Inclusion and exclusion criteria

The scope of this study was limited to simulations of metalloproteins, specifically including articles that examined the effects of metals and metal ligands on proteins in molecular dynamics simulations conducted using AMBER, CHARMM, GROMACS, and NAMD software. Only articles published in English between 2009 and 2024 were selected, with exclusion criteria applied to books and reviews.

Screening and selection of records

The selection process followed PRISMA guidelines [11], and duplicate removal was performed using Rayyan [12]. Two reviewers screened articles by examining the titles and abstracts (P. P. P. and N. F.), then those selected were assessed by a third reviewer (A. A.) through discussion.

Data extraction and presentation

Metalloprotein simulations were examined based on the applied parameterization methods, including molecular mechanics calculations, semi-empirical methods, and Density Functional Theory (DFT). Additionally, data obtained were organized into tables, and each literature was identified by author name and publication year.

Quality appraisal

The criteria for quality appraisal presented and explained in table 1 were used to determine article inclusion by selecting 'Yes' or 'No' indicators to show the adherence of each study. A team of reviewers, consisting of P. P. P, N. F., A. D, and A. A, conducted the quality assessment of articles.

Table 1: Criteria appraised to evaluate the quality of the study that was included

Criteria Description
Study design rationale The study design aligns with the focus on modeling atomic interactions featuring metals in chemistry and biochemistry.
Reproducibility Methods for conducting molecular dynamics simulations, including parameter utilization, are clearly described. Materials and samples are also presented in detail.
Replication Experiments are conducted using molecular dynamics simulations in triplicate or more to ensure reliability.
Negative/positive control Results are compared with existing models or experimental data to validate the accuracy of the simulations.
Parameterization of metal proteins Investigation methods for parametrizing metals in proteins, regardless of the simulation context.
Study adequacy The study design is sufficient to show the mechanisms underlying metal interactions in molecular dynamics simulations.

RESULTS AND DISCUSSION

The screening process was performed using the PRISMA method to select 550 articles, among which 60 were included for review. These 60 comprised full studies focusing on the analysis of metal parameters in molecular dynamics simulations. Specifically, the studies were limited to four software platforms, including Bonded Model Metal Center Parameter Builder, non-bonded metal model, and two types of molecular dynamics simulations. The systematic review method of this study has some limitations concerning the literature search bias, as it depended only on Science Direct and PubMed. In addition, both inclusion and exclusion criteria could lead to biasing because they based themselves on the diversity backgrounds of the four authors albeit within the pharmaceutical field.

The parameter builder for metal centers in bonded model

A Bonded Model in molecular dynamics simulations refers to metals bound to ligands in proteins. Preparation and parameterization of bonded model with QM/MM are necessary before conducting molecular dynamics simulations. Different software packages, such as MCPB, have been developed for the preparation of bonded model. Specifically, MCPB. py uses the Python programming language to generate parameters for over 80 metal ions and force fields, enabling the replication of ligands and metal complexes [13].

The first step in preparing simulations with MCPB. py is creating an input file for parameterization using DFT with the B3LYP/6-31G* basis set. The input file is compatible with Gaussian and GAMESS-US software, and this initial step produces model data cut into small, standard, and large fingerprints. The resulting input file passes through three stages including geometry optimization, Force Constant (FCC) calculation, and Controlled Electrostatic Potential (RESP) calculation, which respectively describe the ideal geometry, bond strength, and charge distribution. The second step comprises using the Seminario method to derive force field parameters. Alternatively, Z-matrix and empirical methods are available, providing flexibility for the parameter derivation process. The empirical method does not require Gaussian or GAMESS-US calculations for force constant determination but relies on Gaussian calculations for RESP values and only supports zinc ion modeling. In the third step, RESP charges are computed and matched with ligand and complex metal residue data. Charging options include calculating charges for both metals and ligands in a single residue unit or solely for selected backbone atoms using specific force fields. In the fourth step, a new Protein Data Bank and input files are obtained along with force fields for proteins (AMBER FF19SB) and ligands (GAFF2). Additionally, ionization data and solvation model for molecular dynamics simulations are generated. To perform further simulations, MCPB. py can be seamlessly integrated with software such as AMBER [14], CHARMM [15], GROMACS [16], and NAMD [17].

The application of MCPB for parameterizing Atovaquone and PfCytb-ISP proteins included validating the Fe2+ parameters, which described b-type hemes and the [2FE-2S] cluster. These parameters were verified to effectively maintain the metal ions in the coordination sphere during all-atom MD simulations of PfCytb-ISP protein present in a phospholipid bilayer. The validated parameters proved beneficial not only for PfCytb-ISP protein but also for other metalloproteins sharing similar coordination environments [18]. Additionally, the Restrained Electrostatic Potential (RESP) charges for the metal ion, sulfur atoms, as well as interacting atoms were computed to complete the parameter and coordinate files for cofactors and all associated residues.

Table 2: Systematically retrieved studies focusing on molecular dynamics simulations of metalloproteins using software such as AMBER, CHARMM, GROMACS, and NAMD

Ref. The simulated metal proteins The metals being analyzed Parameterization of metals Simulation software
[23] The coordination structures of metal ions in the NMR-derived structures. Zn2+ MCPB was used for metal parameterization. Force field parameters for Zn(II) were derived from DFT calculations at the TPSSKCIS/Zn=LANL2DZ/6-31G* level of theory. AMBER
[25] The influence of putidaredoxin (Pdx) on the conformation of the P450cam-CN complex Fe3+ Fe metal parameterization was performed with MCPB. The molecular electrostatic potential (MEP) was calculated using the HF/6-31G* level of theory. AMBER
[26] Virtual screening for the Indoleamine 2,3-dioxygenase 1 (IDO1) target, simulation using classical molecular dynamics Fe Parameters were obtained from quantum mechanics computations using B3LYP, and restrained electrostatic potential (RESP) charges were determined through fitting with the RESP method using ChgModB. MCPB is applied for metal parameterization. AMBER
[27] Electron Transfer to Heme in Cytochrome P450, combined with Umbrella Sampling Fe2+ Using MCPB for metal parameterization. Partial atomic charges were acquired through the RESP method at the HF/6-31G* theory level. AMBER
[28] The influence of ligands on apyrimidinic endonuclease 1 Mg2+, Ni2+ Using B3LYP-D3/6-311G (d, p) to construct convergent side-chain models and MCPB for metal parameterization. AMBER
[21] Molecular dynamics of Metallo-β-lactamase with drugs such as ampicillin, imipenem, nitrocefin, and aztreonam Zn2+ Application of the software GAMESS to obtain metal coordination centers and electrostatics. Additionally, PM3 and DFTB3 are used to determine the relaxed bonded model and the metal is parametrized with MCPB. AMBER
[29] AMBER force field for parameterization of a set of cobalt-containing Systems for bonded and non-bonded model Co B3LYP/6-311++G** level of theory. AMBER
[30] The electron-transporting iron-sulfur protein rubredoxin, combined with Quantum mechanics-discrete molecular dynamics QM/DMD Mg2+ DFT was used for the QM region description, while the BP86 method enhanced with advanced features such as identityresolution, empirical dispersion correction, and a double zeta quality basis set, was applied for single-point energy calculations. COSMO continuum solvation model was used todetermine solvent effects, with a dielectric constant. All QM calculations were conducted with TURBOMOLE v.6.3. DMD, a modified type of MD, applied an implicit solvent atomistic model of proteins based on CHARMM/EEF1 FF. AMBER
[31] The use of Gaussian accelerated molecular dynamics (GaMD) for iron-binding proteins in Botrytis cinerea Fe3+ The calculation of binding energy is performed using DFT. AMBER
[32] Homologue of an Extradiol-cleaving Catecholic Dioxygenase, combined with Stochastic boundary molecular dynamics (SBMD) Fe2+ The substrate DHBP force-field parameters were determined through SwissParam, with additional optimization of partial charge parameters performed using Gaussian09. For the iron and dioxygen groups, force-field parameterizations were manually added to the topology file, while hydrogen atoms for the BphC enzyme were incorporated using the CHARMM36 force field through the HBUILD module of the CHARMM package. CHARMM
[33] Metallothioneins simulation using Steered molecular dynamics (SMD) Zn+2 The Zn7MT2 system was generated by substituting four Cd(II) ions with Zn(II). PROPKA was used to assign protonation states of side chains at pH 7.0, except for Cys residues, which were deprotonated. Protein and Cys residues were modeled using the AMBER FF19SB force field, alongside recently published cysteine-Zn(II) force field parameters. GROMACS
[34] Ferritins from Gram-Negative Bacteria simulation using Random acceleration molecular dynamics (RAMD) Fe2+ The CHARMM36 Force Field wasused to define parameters for non-bonded ferrous and ferric ions. NAMD

Parameterization was conducted for the Cytochrome P450cam-CN in Complex with Putidaredoxin to examine the iron-sulfur cluster of Pdx and the heme iron moiety for ferric camphor-bound, camphor-free and cyanide-bound heme. The focus was on analyzing Distance Measurements of P450cam in different states between residues 48 and 179 in the crystal structures of reduced Cytochrome P450cam-putidaredoxin complex (PDB ID 4jx1), substrate-free P450cam at low [K+] (PDB ID 3L62), and camphor-bound P450cam at low [K+] (PDB ID 3L63) [19]. MCPB can be used to parameterize myoglobin, which contains Fe2+ metal ions bound to histidine [20]. Virtual screening method for heme-containing proteins could be combined with MCPB for parameterization, using B3LYP functional calculations and Grimme D3 dispersion correction alongside Becke-Johnson (BJ)-damping. Additionally, charges of the backbone heavy atoms (ChgModB) were incorporated in the RESP calculations and decomposition analysis of Wiberg bond order in Natural Atomic Orbital (NAO) basis for Fe-O bond was performed to evaluate total bond formation [21]. The functional pathway was evaluated to provide insight into the hydrogen bonding features that define proton transport pathways in [FeFe]-hydrogenase [22].

Combining bonded and non-bonded model simulations is a viable method, requiring careful consideration of the electrostatic bonding strength between zinc ions and deprotonated amino acids such as cysteine and aspartate to stabilize the metal ion at the active site of a protein. The use of MCPB. py in the limited model maintains the coordinated positions of the zinc ions, leading to simulations similar to the crystallographic results [21].

A precise representation of the allosterically inhibited zinc-bound state required fine-tuning zinc-ion coordination in S. aureus CzrA protein [23]. The effect of Putidaredoxin on the Cytochrome P450cam-CN complex, with palladium as the metal center in the ligand, was studied to understand the Pdx-induced conformational changes that promote compound I formation and product release [24]. Moreover, the use of MCPB can determine electrostatic potential and charge distribution, aiding metal parameterization in protein simulations. Observations of structural changes obtained through molecular dynamics simulations provide insight into the conformation of systems featuring specific molecules such as camphor, the P450cam enzyme, and the Pdx protein [24]. Furthermore, three types of conformations were identified, namely closed, intermediate, and open. Pdx binding to the CN−/camphor complex triggered major changes in the intermediate conformation, including peptide bond reversal at Gly248 and the effect on Thr252 reorientation, which affected the width of the helix I protrusion. In the intermediate conformation, it also allowed the formation of hydrogen bonds with CN−, Thr252, and Val247, appearing as catalytic water in the ferrous−O2 complex crystal structure [25].

Virtual screening helps to identify novel inhibitors of Indoleamine 2,3-dioxygenase 1 (IDO1), an intracellular heme-containing dioxygenase associated with immunosuppressive effects in cancer. In this context, MCPB was used to parameterize the heme site, facilitating modeling and integration with classical force fields. Additionally, quantum mechanics (QM) was performed using Gaussian software [26], followed by computing atomic partial charges on heme states and conducting RESP fitting with ChgModB [35]. Electrostatic potential analysis aids in understanding the influence of metals on ligands, with heme parameters available from the Bryce group database [36] and the RESP ESP charge database (R. E. DD. B) [37].

The construction of the Mg-O force field was accomplished with MCPB due to the highly fluctuating nature of the Mg2+ion. The trajectory of Mg2+ was observed through snapshotting at the equilibrium point, while metal active site movement in the system was visualized by monitoring bond lengths during simulations [28]. Combining MD simulations, umbrella sampling simulations, and electron transfer pathways analysis, the mechanisms of long-range electron transfer between the [2Fe-2S]-containing redox partner and the heme domain in full-length P450TT were deciphered [27]. Metalloproteins often experience energetically unfavorable conformational changes or binding events over extended timescales. Umbrella sampling addresses these challenges by guiding the simulation towards specific states of interest.

Parameterization of metals using a non-bonded model

A non-bonded model generally refers to a metal not bound to a ligand but present in a protein. Moreover, this review identified specific AMBER force field parameters available for 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 used to determine cobalt bond stretching and angle bending parameters for the non-bonded model. Additionally, the charge model consists of RESP and formal charges [38]. The movement of a non-bonded model can be monitored through RMSD values during vacuum and water simulations, allowing for observations of changes in cobalt, such as bond length, bond strain constant, and bond equilibrium [29]. Free Energy Landscape (FEL) analysis can be used to see metal ion mechanisms on conformational distribution and complex stability [39].

Preparation of simulations for cadmium includes using the 12-6 Lennard-Jones (LJ) non-bonded model, particularly suitable for divalent cations, with water model such as TIP3P, SPC/E, and TIP4PEW [40, 41]. The non-bonded model treats the metal ion as a particle interacting non-covalently with surrounding molecules. Coordination bonds in this model are governed by electrostatic and van der Waals potentials. Different non-bonded models apply distinct strategies to handle electrostatic interactions, ensuring accurate coordination geometry description [42, 9].

Proteins and metals are prepared after parameterization using the tleap program and the AMBER99SB force field, then simulated with AMBER software [43, 14]. Monovalent metal simulations are enhanced with optimized parameters for the 12–6 lJ non-bonded model [44] and LJ 12-6 Potential parameterization [45]. Meanwhile, divalent model parameters can be analyzed to determine effective ion radii, diffusion constants, water exchange rates, and ion–water interactions [46].

The 12-6-4 lJ Type Non-bonded Model is suitable for simulations of divalent [47], trivalent, and tetravalent metals [48]. The binding of Cm(III) and Th(IV) with Human Transferrin, representing non-bonded types, can be parameterized by examining the dynamic hydration shell structure in an aqueous solution [49]. Additionally, umbrella sampling is used to analyze metal-protein systems, enabling the evaluation of the Potential mean Force (PMF) of ion release [50]. The parameterization of Highly Charged Metal Ions with the 12-6-4 lJ model can be conducted based on a previous study [1]. Panteva et al. recently found that the combination of 12-6-4 with the SPC/E water model produced the best performance for Mg(II) among 17 different non-bonded models investigated [51].

Comparison was made between zinc parameters in the AMBER force field (ZAFF) and a recently developed non-bonded force field (NBFF) to assess the accuracy in reproducing zinc(II)-protein dynamics. To facilitate this process, six zinc fingers were selected as benchmarks due to the architectural heterogeneity, binding mode, function, and reactivity. Additionally, order parameters (S2) of all backbone N-H bond vectors were computed for each system through repeated MD simulations [52]. Geometric and electronic behavior of the midodrine complex with iron(III) using DFT with the B3LYP/6–311G++(d,p) method for the ligand and B3LYP/lANL2D for this metal was carried out to see the molecular orbital and electrostatic potential [53]. To understand the nature of the cofactor, it is very important to analyze the Van der Waals interactions in bridging sulfur and iron surfaces [54]. Apart from the DFT method, semi-empirical methods such as ZINDO can also be used to parameterize metalproteins because they are very accurate in describing different spin and electronic spectra [55].

Types of molecular dynamics simulations conducted in modeling metalloproteins

Force field parameters were established for the active site of LOX containing Fe2+, representing the bonded model. The construction of the active site applied a quantum mechanical (QM) method with geometry optimization using B3LYP/6-31G(d,p) due to the ability to produce equilibrium bond lengths and angles between Fe(II) ion and amino acids. Additionally, force constants were computed by approximating the Hessian matrix. Charges were calculated using electrostatic surface potential (ESP) with the Merz-Kollman Scheme (MK) method, while electrostatic interactions were handled with the restrained electrostatic surface potential (RESP). Following the parameterization process, molecular dynamics simulations were conducted using AMBER software, during which bond distances between Fe and amino acids, as well as water molecules, were monitored [56].

The QM/MM method was implemented for Protoporphyrin IX in Human Ferrochelatase using ComQum software, a program applying Turbomole to build a system consisting of the protoporphyrin ring, ferrous iron, coordinating water molecules, and relevant residues treated at the QM level [57-59]. Subsequently, minimization steps were performed with the QM region computing wave functions and polarization using point charges from the AMBER library. Protein preparation was conducted with AMBER and protonation; then Quantum calculations were performed using DFT with the DZpdf basis set. Fe was computed using the 6-31G(d) basis set, while more accurate single-point energies were estimated through Becke B3LYP functional method and the 6-311+G(2d,2p) basis set. quantum mechanical thermodynamic cycle perturbation (QTCP) proved effective in computing proton transfer at the metal site of macromolecules [60]. QM/MM modeling was used to observe the reaction coordinates from the intermediate state to the Fe2+ transition down to the electron level. Simulations were continued to analyze the protein and metal atomic levels with a Molecular Mechanics (MM) method, while QM/MM energies were collected from snapshot trajectories [61]. Additionally, Fe metal applied the LANL2DZ basis set, and the QM region was constructed using the hybrid functional B3LYP with the 6-31G(d) basis set for the 5-lipoxygenase region [62]. Electronic parameter calculations were performed following the MNDO99 method for semi-empirical calculations. Treatment on proteins used the CHARMM22 and CHARMM27 force fields, while adjustments were made with the DL_POLY module in Chemshell. Electronic parameter computations were analyzed with QM/MM to investigate hydrogen bond charges and connectivity. Considering the transition state was found to be crucial in modeling metalloproteins [63].

Applying the metal molecular fractionation with the conjugate caps (metal-MFCC) method, metalloproteins can be accurately modeled by using M06-2X/6-31G(d) level quantum system computations with linear potential energy and atomic force calculations [64]. Furthermore, the total energy of the metalloprotein is determined by adding the interaction energies between nonadjacent residues in close spatial contact, the energies of metal-binding groups, and the QM energies of neighboring residues, which is a highly suitable method for Zn2+, Cu2+, and Cu+ [26]. To evaluate the metal effects on the active site, average distances between Zn(II) and the donor atoms along the QM/MM MD simulation are used [65].

The data obtained from the use of the Broyden−Fletcher−Goldfarb−Shanno (L-BFGS) algorithm to optimize the geometry in an energy minimization search are essential for observing the distance between the reactants of the triplet, quintet, and septet states of the Catalytic Mechanism for the 2,3-Dihydroxybiphenyl ring. Since the QM region coincides with the electron density distribution, it is essential to study it, particularly when examining Fe [32, 66]. Furthermore, the CHARMM-METAL force field can be used to analyze the properties of Fe metal [67].

The multiscale method known as our own N-layered integrated molecular orbital and molecular mechanics (ONIOM) integrating DFT with molecular mechanics is used to study Cu and Zn metals in metalloproteins [68]. Additionally, a semi-empirical method referred to as GFN2-xTB can be used to investigate effectiveness and improve local structure dependability in metalloproteins [69]. Through ONIOM, electronic properties and interactions in metalloproteins can be more accurately described by merging multiple levels of theory, such as DFT and semi-empirical [70]. The combination of DFT simulations with geometry optimization and energy calculations forms an additional method for observing the transition state and examining the catalytic process of L-tyrosine and Cu2+ [71].

Torsional parameter dependence on partial charges is a classic issue with histidine-metal interactions, mainly including Zn and Cu. DFT calculations were performed to address this issue by using a triple-ζ polarized basis set (def2-TZVP) for metals and a split-valence polarized basis set (def2-SVP) for nonmetal atoms [72]. Furthermore, point charges were calculated with TURBOMOLE using the RESP method [59] and ChemShell [73]. Moreover, the significance of parameterizing torsional terms was shown by molecular dynamics simulations, specifically when focused on the metal site or performing free-energy calculations.

The iron dioxygen complex parameters are derived using the heme parameters found in AMBER ff9X and AMBER ff14SB force fields [74, 75]. The ORCA Quantum Chemistry Program was used to optimize the QM aspects of the model systems. However, determining the Fe-S parameters provided specific challenges in molecular dynamics simulations, which included the use of RESP for frequency and charge calculations as well as B3LYP for geometry optimization. The Seminario method was then applied to identify parameters associated with group or atom interactions [76]. Parametrization of Iron−Sulfur interactions was conducted through DFT with the TZP basis set, followed by point cargo calculations in the gas phase using HF/631G* and RESP calculations [77].

The influence of force fields (FF) on metalloproteins, such as CHARMM36, varied in the ability to describe interactions with monovalent ions. However, the Drude FF and DFTB3/3OB method showed superior performance, specifically for Ca2+ binding sites, signifying the importance of considering polarization effects for divalent cations [78]. QM/MM calculations could be conducted post-molecular dynamics by extracting Pb2+ trajectory data from connected amino acids. Solvent exposure of this Pb metal was evaluated using solvent-accessible surface calculations. The QM region was ass through DFT with the LACVP** basis set, incorporating an effective core potential for Pb2+, and atomic charges were calculated using B3LYP/lACV3P**++level [79]. This analysis provided information on the distance between the metal and protein, along with the resulting energy in the gas phase.

Random acceleration molecular dynamics (RAMD) simulations were used to analyze ferrous iron and dioxide interactions in metalloproteins, providing a cost-effective exploration of phase space [34]. The atomic multipole-optimized energetics for biomolecular simulation (AMOEBA) force field enhanced simulation accuracy for Mg2+/Ca2+ in a variety of metalloproteins, particularly in systems comprising highly charged molecules [80]. The study of prolyl oligopeptidase POP-Rh2 cyclopropanase included DFT calculations combined with Grimme D3 Correction for the solvent dispersion model. Molecular docking using the GaudiMM algorithm played a crucial role in examining the coordination of the rhodium cofactor with amino acid donors [81]. Subsequent molecular dynamics simulations clarified carbene orientation and binding route conformations, further analyzed with GpathFinder [82]. The assessment of styrene local interaction was conducted through molecular docking using GaudiMM and Gold [83]. However, due to the complexity of the study, it is importance of comprehensive modelling that includes long-term mutation effects [84].

Force fields could be tailored for specific metal coordination geometries and ligand interactions. Quantum mechanics-discrete molecular dynamics (QM/DMD) simulations offered a balanced description of metalloprotein structure, dynamics, and electronic structure for rubredoxin [30]. Moreover, transition metal-containing enzyme mechanisms, particularly in redox steps, could be explained with a saturated basis set. DFT was used in the description of the QM region due to being a robust and reliable method for simulating biologically relevant metals [71, 85]. Single-point calculations were performed with various basis sets using Jaguar [86] and Gaussian programs [85].

The poisson−boltzmann (PB) continuum solvent method, used for simulating the aqueous medium [87], permits the extraction of data from molecular dynamics simulations at specific frames. In these simulations, QM calculations were conducted for each fragment, with other residues considered as background charges and solvation effects incorporated by solving the PB equation using the Delphi program [88]. Hybrid QM/MM calculations on metal-bound structures, such as gol B, facilitate the characterization of the most preferred protonation state for gold binding motifs and the assessment of structural features influencing Au(I) coordination in the protein. The QM region, treated at the DFT level with the LACVP** basis set and the B3LYP functional, provides a comprehensive description of geometries and reaction profiles for transition metal-containing compounds, including gold(I). Meanwhile, the MM system, described with the OPLSA_2001 force field, consists of protein atoms and a 10 Å shell of water molecules [87].

Single-point energy calculations use a larger basis set, namely 6–311+G(2d,p), for non-silver atoms and an Effective Core Potential (ECP), such as Def2-TZVPD, for silver atoms. Additionally, ECPs replace core electrons with a pseudopotential, enhancing computational efficiency in metal simulations [89, 90]. The atom-bond electronegativity equalization method (ABEEMσπ) fluctuating charge polarizable force fields provide a medium for parameterizing charge distributions in metalloproteins containing transition metal atoms. This also includes reference charges, valence state electronegativities, and valence state hardnesses [91].

Analysis of Silver-Lactoferrin Nanocomplexes requires molecular dynamics simulations to observe amino acids directly or indirectly interacting with silver while identified binding sites correspond to amino acids that often bind metal ions. DFT is used to evaluate the mechanism of silver ion reduction to elemental silver, followed by an examination of the localization of highest occupied (HOMO) and lowest unoccupied (LUMO) molecular orbitals [92].

The role of Metal Cofactors in azurin is assessed through temperature replica exchange molecular dynamics (REMD) on both apo and holo forms. REMD facilitates the observation of unfolding processes and transition temperatures, providing insights into conformational ensembles and unfolding pathways [93]. Combining cysteine-Zn(II) parameters from the AMBER force field [89] with steered molecular dynamics (SMD) simulations allows for the analysis of metal effects on Mammalian Metallothioneins (MTs), clarifying intra-and interdomain interactions as well as the implications for zinc buffering properties. Understanding the system folding mechanism, zinc site stabilities, and cellular zinc buffer components is crucial for comprehensive analysis [33]. The use of SMD can analyze how the presence of metallic ligands alters the unfolding dynamics by forming strong covalent bonds with cysteine sulfur atoms [94]. Gaussian accelerated molecular dynamics (GaMD) is an advanced computational sampling method that enhances efficiency by introducing harmonic potentials [95]. This reduces energy barriers and smooths out the system's potential energy surface. Following this, binding energy calculations are performed by randomly selecting five structures from a 150 ns trajectory and using DFT for these calculations [31].

CONCLUSION

In conclusion, this study found that modeling and parameterization of metal proteins were categorized into Bonded and Non-Bonded Model. The Bonded Model incorporates MCPB, a Python-based software that facilitates parameter construction for over 80 metal ions and force fields in molecular dynamics simulations. MCPB uses DFT to generate input files for parameterization, derives force field parameters through the Seminario method, and matches RESP charges with ligand and complex metal residue data. This can be integrated with software, including AMBER for further simulations, as well as applied for virtual screening and understanding metal-ligand interactions. Meanwhile, the non-bonded model evaluates metals in proteins, such as cobalt and cadmium, through AMBER force field parameters and the Seminario method. The 12-6 lJ non-bonded model is suitable for divalent, trivalent, and tetravalent metals, with Zinc parameters being compared for accuracy. Additionally, the force fields suitable for modeling unbound metal proteins include AMBER FF19SB, FF14SB, and ff9X, along with CHARMM36, CHARMM22, and CHARMM27, as well as CHARMM-Metal. Molecular dynamics simulations can be conducted using various methods, such as classical molecular dynamics, Umbrella Sampling, QM/DMD, SBMD, SMD, GaMD and RAMD.

ACKNOWLEDGMENT

The authors are grateful to LPPM Universitas Andalas for the support, guidance, and financial assistance rendered.

FUNDING

Generous funding was received from LPPM Universitas Andalas under contract number T/137/UN.16.17/PT.01.03/KO-RPT/2022.

AUTHORS CONTRIBUTIONS

P. P. P. conceptualized this study, then conducted methodology development, validation, formal analysis, original draft preparation, review, and editing alongside two other authors, including N. F., and A. A. Additionally, software utilization and investigation were performed by P. P. P, N. F., A. A., and A. D. Visualization and supervision were handled by P. P. P and N. F., while funding assistance was acquired by P. P. P, N. F., and A. D.

CONFLICT OF INTERESTS

The authors declare no conflicts of interest.

REFERENCES

  1. Li P, Song LF, Merz KM. Parameterization of highly charged metal ions using the 12-6-4 lJ-type nonbonded model in explicit water. J Phys Chem B. 2015;119(3):883-95. doi: 10.1021/jp505875v, PMID 25145273.

  2. Zhang Y, Zheng J. Bioinformatics of metalloproteins and metalloproteomes. Molecules. 2020;25(15):1-23. doi: 10.3390/molecules25153366, PMID 32722260.

  3. Chalkley MJ, Mann SI, DeGrado WF. De novo metalloprotein design. Nat Rev Chem. 2022;6(1):31-50. doi: 10.1038/s41570-021-00339-5, PMID 35811759.

  4. Vidossich P, Magistrato A. QM/MM molecular dynamics studies of metal-binding proteins. Biomolecules. 2014;4(3):616-45. doi: 10.3390/biom4030616, PMID 25006697.

  5. Raut K, Kothawade S, Pande V, Bole S, Netane S, Autade K. Design of potent anticancer molecules comprising pyrazolyl thiazolinone analogs using molecular modelling studies for pharmacophore optimization. Asian J Pharm Clin Res. 2023;16(8):84-93. doi: 10.22159/ajpcr.2023.v16i8.47665.

  6. Dinnimath BM, Gowda P, Naik A. Development of organometallic compounds of schiff bases with diverse applications. Int J Pharm Pharm Sci. 2023;15(6):1-15. doi: 10.22159/ijpps.2023v15i6.47362.

  7. Di Felice R, Mayes ML, Richard RM, Williams Young DB, Chan GK, de Jong WA. A perspective on sustainable computational chemistry software development and integration. J Chem Theory Comput. 2023;19(20):7056-76. doi: 10.1021/acs.jctc.3c00419, PMID 37769271.

  8. Rodzik A, Pomastowski P, Sagandykova GN, Buszewski B. Interactions of whey proteins with metal ions. Int J Mol Sci. 2020;21(6):2156. doi: 10.3390/ijms21062156, PMID 32245108.

  9. Li P, Merz KM. Metal ion modeling using classical mechanics. Chem Rev. 2017;117(3):1564-6. doi: 10.1021/acs.chemrev.6b00440, PMID 28045509.

  10. Iqhrammullah M, Rizki DR, Purnama A, Duta TF, Harapan H, Idroes R. Antiviral molecular targets of essential oils against SARS-CoV-2: a systematic review. Sci Pharm. 2023;91(1):15. doi: 10.3390/scipharm91010015.

  11. Page MJ, McKenzie JE, Bossuyt PM, Boutron I, Hoffmann TC, Mulrow CD. The PRISMA 2020 statement: an updated guideline for reporting systematic reviews. BMJ. 2021;372:71. doi: 10.1136/bmj.n71.

  12. Ouzzani M, Hammady H, Fedorowicz Z, Elmagarmid A. Rayyan a web and mobile app for systematic reviews. Syst Rev. 2016;5(1):210. doi: 10.1186/s13643-016-0384-4, PMID 27919275.

  13. Li P, Merz KM. MCPB. py: a python-based metal center parameter builder. J Chem Inf Model. 2016;56(4):599-604. doi: 10.1021/acs.jcim.5b00674. PMID 26913476.

  14. Case DA, Aktulga HM, Belfon K IY, Ben Shalom JT, Berryman SR, Brozell DS, Cerutti TE. Izad and PAK. Amber 2023. San Fr, University California; 2023.

  15. Brooks BR, Brooks CL, Mackerell AD, Nilsson L, Petrella RJ, Roux B. CHARMM: the biomolecular simulation program. J Comput Chem. 2009;30(10):1545-614. doi: 10.1002/jcc.21287, PMID 19444816.

  16. Abraham MJ, Murtola T, Schulz R, Pall S, Smith JC, Hess B. Gromacs: high-performance molecular simulations through multi-level parallelism from laptops to supercomputers. Software X. 2015;1-2:19-25. doi: 10.1016/j.softx.2015.06.001.

  17. Phillips JC, Hardy DJ, Maia JD, Stone JE, Ribeiro JV, Bernardi RC. Scalable molecular dynamics on CPU and GPU architectures with NAMD. J Chem Phys. 2020;153(4):044130. doi: 10.1063/5.0014475, PMID 32752662.

  18. Chebon Bore L, Sanyanga TA, Manyumwa CV, Khairallah A, Tastan Bishop OT. Decoding the molecular effects of atovaquone-linked resistant mutations on plasmodium falciparum cytb-isp complex in the phospholipid bilayer membrane. Int J Mol Sci. 2021;22(4):1-32. doi: 10.3390/ijms22042138, PMID 33670016.

  19. Lee YT, Wilson RF, Rupniewski I, Goodin DB. P450cam visits an open conformation in the absence of substrate. Biochemistry. 2010;49(16):3412-9. doi: 10.1021/bi100183g, PMID 20297780.

  20. Aghaie Kheyrabadi F, Shareghi B, Farhadian S. Molecular interaction between cypermethrin and myoglobin: spectroscopy and molecular dynamics simulation analysis. J Mol Liq. 2024;395. doi: 10.1016/j.molliq.2024.123949.

  21. Eshtiwi AA, Rathbone DL. A modified bonded model approach for molecular dynamics simulations of New Delhi metallo-β-lactamase. J Mol Graph Model. 2023;121:108431. doi: 10.1016/j.jmgm.2023.108431, PMID 36827734.

  22. Ginovska Pangovska B, Ho MH, Linehan JC, Cheng Y, Dupuis M, Raugei S. Molecular dynamics study of the proposed proton transport pathways in [FeFe]-hydrogenase. Biochim Biophys Acta. 2014;1837(1):131-8. doi: 10.1016/j.bbabio.2013.08.004, PMID 23981729.

  23. Chakravorty DK, Wang B, Lee CW, Guerra AJ, Giedroc DP, Merz KM. Solution NMR refinement of a metal ion bound protein using metal ion inclusive restrained molecular dynamics methods. J Biomol NMR. 2013;56(2):125-37. doi: 10.1007/s10858-013-9729-7, PMID 23609042.

  24. Hollingsworth SA, Poulos TL. Molecular dynamics of the P450cam-Pdx complex reveals complex stability and novel interface contacts. Protein Sci. 2015;24(1):49-57. doi: 10.1002/pro.2583, PMID 25307478.

  25. Chuo SW, Wang LP, Britt RD, Goodin DB. An intermediate conformational state of cytochrome P450cam-CN in complex with putidaredoxin. Biochemistry. 2019;58(18):2353-61. doi: 10.1021/acs.biochem.9b00192, PMID 30994334.

  26. Zou Y, Hu Y, Ge S, Zheng Y, Li Y, Liu W. Effective virtual screening strategy toward heme-containing proteins: identification of novel IDO1 inhibitors. Eur J Med Chem. 2019;184:111750. doi: 10.1016/j.ejmech.2019.111750, PMID 31610376.

  27. Wang Z, Shaik S, Wang B. Conformational motion of ferredoxin enables efficient electron transfer to heme in the full length P450TT. J Am Chem Soc. 2021;143(2):1005-16. doi: 10.1021/jacs.0c11279, PMID 33426875.

  28. Wu Z, Duan H, Cheng Y, Guo D, Peng L, Hu Y. A novel ligand swing mediated active site coordination change of human apurinic apyrimidinic endonuclease 1: a potential cytotoxic mechanism of nickel ion in the base excision repair. Chem Phys. 2022;555. doi: 10.1016/j.chemphys.2022.111456.

  29. Haghshenas H, Tavakol H, Kaviani B, Mohammadnezhad G. AMBER force field parameters for cobalt containing biological systems: a systematic derivation study. J Phys Chem B. 2020;124(5):777-87. doi: 10.1021/acs.jpcb.9b10739, PMID 31912730.

  30. Sparta M, Shirvanyants D, Ding F, Dokholyan NV, Alexandrova AN. Hybrid dynamics simulation engine for metalloproteins. Biophys J. 2012;103(4):767-76. doi: 10.1016/j.bpj.2012.06.024, PMID 22947938.

  31. Wu R, Xie D, Du J. The binding pattern of ferric iron and iron binding protein in botrytis cinerea. Comput Biol Med. 2024;178:108686. doi: 10.1016/j.compbiomed.2024.108686, PMID 38850956.

  32. Wang J, Chen J, Tang X, Li Y, Zhang R, Zhu L. Catalytic mechanism for 2,3-dihydroxy biphenyl ring cleavage by nonheme extradiol dioxygenases BphC: insights from QM/MM analysis. J Phys Chem B. 2019;123(10):2244-53. doi: 10.1021/acs.jpcb.8b11008, PMID 30776233.

  33. Singh AK, Pomorski A, Wu S, Peris Diaz MD, Czepczynska Krezel H, Krezel A. The connection of α- and β-domains in mammalian metallothionein-2 differentiates Zn(II) binding affinities affects folding and determines zinc buffering properties. Metallomics. 2023;15(6):29. doi: 10.1093/mtomcs/mfad029, PMID 37147085.

  34. Pietra F. New vistas on the recruiting of ferrous iron and dioxygen by ferritins: a case study of the escherichia coli 24-mer ferritin by all-atom molecular dynamics in aqueous medium. Chem Biodivers. 2018;15(8):e1800197. doi: 10.1002/cbdv.201800197, PMID 29873188.

  35. Peters MB, Yang Y, Wang B, Fusti Molnar L, Weaver MN, Merz KM. Structural survey of zinc-containing proteins and the development of the zinc amber force field (zaff). J Chem Theory Comput. 2010;6(9):2935-47. doi: 10.1021/ct1002626, PMID 20856692.

  36. AMBER parameter database. In: Manchester UK: University of Manchester. Available from: http://research.bmh.manchester.ac.uk/bryce/amber. [Last accessed on 14 Aug 2024]

  37. Dupradeau FY, Cezard C, Lelong R, Stanislawiak E, Pecher J, Delepine JC. A database for RESP and ESP atomic charges and force field libraries. Nucleic Acids Res. 2008;36Suppl 1.

  38. Donati G, D’Amore VM, Russomanno P, Cerofolini L, Amato J, Marzano S. Theoretical and experimental studies on the interaction of biphenyl ligands with human and murine PD-L1: up-to-date clues for drug design. Comput Struct Biotechnol J. 2023;21:3355-68. doi: 10.1016/j.csbj.2023.06.006, PMID 37384351.

  39. He Y, Jiang Z, Zeng M, Cao S, Wu N, Liu X. Unraveling potential mechanism of different metal ions effect on anammox through big data analysis molecular docking and molecular dynamics simulation. J Environ Manage. 2024;352:120092. doi: 10.1016/j.jenvman.2024.120092, PMID 38232596.

  40. Li P, Roberts BP, Chakravorty DK, Merz KM. Rational design of particle mesh ewald compatible Lennard Jones parameters for +2 metal cations in explicit solvent. J Chem Theory Comput. 2013;9(6):2733-48. doi: 10.1021/ct400146w, PMID 23914143.

  41. Singh G, Tripathi S, Shanker K, Sharma A. Cadmium-induced conformational changes in type 2 metallothionein of medicinal plant coptis japonica: insights from molecular dynamics studies of apo partially and fully metalated forms. J Biomol Struct Dyn. 2019;37(6):1520-33. doi: 10.1080/07391102.2018.1461688, PMID 29624115.

  42. Popelier PL. Non-covalent interactions from a quantum chemical topology perspective. J Mol Model. 2022;28(9):276. doi: 10.1007/s00894-022-05188-7, PMID 36006513.

  43. Lindorff Larsen K, Piana S, Palmo K, Maragakis P, Klepeis JL, Dror RO. Improved side chain torsion potentials for the amber ff99SB protein force field. Proteins. 2010;78(8):1950-8. doi: 10.1002/prot.22711, PMID 20408171.

  44. Li P, Song LF, Merz KM. Systematic parameterization of monovalent ions employing the nonbonded model. J Chem Theory Comput. 2015;11(4):1645-57. doi: 10.1021/ct500918t, PMID 26574374.

  45. Qiu Y, Jiang Y, Zhang Y, Zhang H. Rational design of nonbonded point charge models for monovalent ions with Lennard Jones 12-6 potential. J Phys Chem B. 2021;125(49):13502-18. doi: 10.1021/acs.jpcb.1c09103, PMID 34860517.

  46. Zhang Y, Jiang Y, Peng J, Zhang H. Rational design of nonbonded point charge models for divalent metal cations with lennard jones 12-6 potential. J Chem Inf Model. 2021;61(8):4031-44. doi: 10.1021/acs.jcim.1c00580, PMID 34313132.

  47. Li Z, Song LF, Li P, Merz KM. Systematic parametrization of divalent metal ions for the OPC3, OPC, TIP3P-FB, and TIP4P-FB water models. J Chem Theory Comput. 2020;16(7):4429-42. doi: 10.1021/acs.jctc.0c00194, PMID 32510956.

  48. Li Z, Song LF, Li P, Merz KM. Parametrization of trivalent and tetravalent metal ions for the OPC3, OPC, TIP3P-FB, and TIP4P-FB water models. J Chem Theory Comput. 2021;17(4):2342-54. doi: 10.1021/acs.jctc.0c01320, PMID 33793233.

  49. Atta Fynn R, Bylaska EJ, Schenter GK, De Jong WA. Hydration shell structure and dynamics of curium(III) in aqueous solution: first principles and empirical studies. J Phys Chem A. 2011;115(18):4665-77. doi: 10.1021/jp201043f, PMID 21500828.

  50. Mishra L, Sundararajan M, Bandyopadhyay T. MD simulation reveals differential binding of Cm(III) and th(IV) with serum transferrin at acidic pH. Proteins. 2021;89(2):193-206. doi: 10.1002/prot.26006, PMID 32892408.

  51. Panteva MT, Giambaşu GM, York DM. Comparison of structural, thermodynamic kinetic and mass transport properties of Mg(2+) ion models commonly used in biomolecular simulations. J Comput Chem. 2015;36(13):970-82. doi: 10.1002/jcc.23881, PMID 25736394.

  52. Bazayeva M, Giachetti A, Pagliai M, Rosato A. A comparison of bonded and nonbonded zinc(II) force fields with NMR data. Int J Mol Sci. 2023;24(6):5440. doi: 10.3390/ijms24065440, PMID 36982515.

  53. El Sayed DS, Khalil TE, Elbadawy HA. Rational and experimental investigation of antihypotensive midodrine-Fe(III) complex: synthesis spectroscopy DFT biological activity and molecular docking. J Mol Struct. 2024;1311. doi: 10.1016/j.molstruc.2024.138421.

  54. Jenney FE, Wang H, George SJ, Xiong J, Guo Y, Gee LB. Temperature-dependent iron motion in extremophile rubredoxins no need for ‘corresponding states’ Sci Rep. 2024;14(1):12197. doi: 10.1038/s41598-024-62261-2, PMID 38806591.

  55. Hagras MA. Respiratory complex iii: a bioengine with a ligand triggered electron-tunneling gating mechanism. J Phys Chem B. 2024;128(4):990-1000. doi: 10.1021/acs.jpcb.3c07095, PMID 38241470.

  56. Moin ST, Hofer TS, Sattar R, Ul-Haq Z. Molecular dynamics simulation of mammalian 15S-lipoxygenase with AMBER force field. Eur Biophys J. 2011;40(6):715-26. doi: 10.1007/s00249-011-0684-5, PMID 21360129.

  57. Ryde U. The coordination of the catalytic zinc in alcohol dehydrogenase studied by combined quantum chemical and molecular mechanics calculations. J Comput Aided Mol Des. 1996;10(2):153-64. doi: 10.1007/BF00402823, PMID 8741019.

  58. Ryde U, Olsson MH. Structure strain and reorganization energy of blue copper models in the protein. Int J Quantum Chem. 2001;81(5):335-47. doi: 10.1002/1097-461X(2001)81:5<335::AID-QUA1003>3.0.CO;2-Q.

  59. Treutler O, Ahlrichs R. Efficient molecular numerical integration schemes. J Chem Phys. 1995;102(1):346-54. doi: 10.1063/1.469408.

  60. Kaukonen M, Soderhjelm P, Heimdal J, Ryde U. Proton transfer at metal sites in proteins studied by quantum mechanical free energy perturbations. J Chem Theory Comput. 2008 Jun 1;4(6):985-1001. doi: 10.1021/ct700347h, PMID 26621239.

  61. Wu J, Wen S, Zhou Y, Chao H, Shen Y. Human ferrochelatase: insights for the mechanism of ferrous iron approaching protoporphyrin IX by QM/MM and QTCP free energy studies. J Chem Inf Model. 2016;56(12):2421-33. doi: 10.1021/acs.jcim.6b00216, PMID 27801584.

  62. Saura P, Marechal JD, Masgrau L, Lluch JM, Gonzalez Lafont A. Computational insight into the catalytic implication of head/tail-first orientation of arachidonic acid in human 5-lipoxygenase: consequences for the positional specificity of oxygenation. Phys Chem Chem Phys. 2016;18(33):23017-35. doi: 10.1039/c6cp03973a, PMID 27489112.

  63. Suardiaz R, Masgrau L, Lluch JM, Gonzalez Lafont A. Introducing mutations to modify the C13/C9 ratio in linoleic acid oxygenations catalyzed by rabbit 15-lipoxygenase: a QM/MM and MD study. ChemPhysChem. 2014;15(18):4049-54. doi: 10.1002/cphc.201402471, PMID 25186660.

  64. Xu M, He X, Zhu T, Zhang JZ. A fragment quantum mechanical method for metalloproteins. J Chem Theory Comput. 2019;15(2):1430-9. doi: 10.1021/acs.jctc.8b00966, PMID 30620584.

  65. Sgrignani J, Magistrato A, Dal Peraro MD, Vila AJ, Carloni P, Pierattelli R. On the active site of mononuclear B1 metallo β-lactamases: a computational study. J Comput Aided Mol Des. 2012;26(4):425-35. doi: 10.1007/s10822-012-9571-0, PMID 22532071.

  66. Hong G, Pachter R. Inhibition of biocatalysis in [Fe-Fe] hydrogenase by oxygen: molecular dynamics and density functional theory calculations. ACS Chem Biol. 2012;7(7):1268-75. doi: 10.1021/cb3001149, PMID 22563793.

  67. Ruiz MP, Aragones AC, Camarero N, Vilhena JG, Ortega M, Zotti LA. Bioengineering a single protein junction. J Am Chem Soc. 2017;139(43):15337-46. doi: 10.1021/jacs.7b06130, PMID 28981262.

  68. Neese F, Wennmohs F, Hansen A, Becker U. Efficient approximate and parallel hartree fock and hybrid DFT calculations. A ’chain-of-spheres’ algorithm for the hartree fock exchange. Chem Phys. 2009;356(1-3):98-109. doi: 10.1016/j.chemphys.2008.10.036.

  69. Bannwarth C, Ehlert S, Grimme S. GFN2-xTB-an accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependent dispersion contributions. J Chem Theory Comput. 2019;15(3):1652-71. doi: 10.1021/acs.jctc.8b01176, PMID 30741547.

  70. Yan Z, Li X, Chung LW. Multiscale quantum refinement approaches for metalloproteins. J Chem Theory Comput. 2021;17(6):3783-96. doi: 10.1021/acs.jctc.1c00148, PMID 34032440.

  71. Jia L, Wang C, Zhang S, Yuan Z, Lu F, Liu Y. Structure-guided protein engineering to reconstruct the biocatalyst properties for efficient l-3,4-dihydroxyphenylalanine production. Chem Eng J. 2023;469:143894. doi: 10.1016/j.cej.2023.143894.

  72. Mera Adasme R, Sadeghian K, Sundholm D, Ochsenfeld C. Effect of including torsional parameters for histidine metal interactions in classical force fields for metalloproteins. J Phys Chem B. 2014;118(46):13106-11. doi: 10.1021/jp5078906, PMID 25410708.

  73. Lu Y, Sen K, Yong C, Gunn DS, Purton JA, Guan J. Multiscale QM/MM modelling of catalytic systems with chemshell. Phys Chem Chem Phys. 2023;25(33):21816-35. doi: 10.1039/d3cp00648d, PMID 37097706.

  74. Shahrokh K, Orendt A, Yost GS, Cheatham TE. Quantum mechanically derived AMBER-compatible heme parameters for various states of the cytochrome P450 catalytic cycle. J Comput Chem. 2012;33(2):119-33. doi: 10.1002/jcc.21922, PMID 21997754.

  75. Ugur I, Chandrasekhar P. Proton relay network in P450cam formed upon docking of putidaredoxin. Proteins. 2020;88(4):558-72. doi: 10.1002/prot.25835, PMID 31597203.

  76. Carvalho AT, Teixeira AF, Ramos MJ. Parameters for molecular dynamics simulations of iron-sulfur proteins. J Comput Chem. 2013;34(18):1540-8. doi: 10.1002/jcc.23287, PMID 23609049.

  77. Carvalho AT, Swart M. Electronic structure investigation and parametrization of biologically relevant iron-sulfur clusters. J Chem Inf Model. 2014;54(2):613-20. doi: 10.1021/ci400718m, PMID 24460186.

  78. Ngo V, Da Silva MC, Kubillus M, Li H, Roux B, Elstner M. Quantum effects in cation interactions with first and second coordination shell ligands in metalloproteins. J Chem Theory Comput. 2015;11(10):4992-5001. doi: 10.1021/acs.jctc.5b00524, PMID 26574284.

  79. Tolbatov I, Re N, Coletti C, Marrone A. Determinants of the lead(II) affinity in pbrR protein: a computational study. Inorg Chem. 2020;59(1):790-800. doi: 10.1021/acs.inorgchem.9b03059, PMID 31829577.

  80. Litman J, Thiel AC, Schnieders MJ. Scalable indirect free energy method applied to divalent cation-metalloprotein binding. J Chem Theory Comput. 2019;15(8):4602-14. doi: 10.1021/acs.jctc.9b00147, PMID 31268700.

  81. Rodriguez Guerra Pedregal J, Sciortino G, Guasp J, Municoy M, Marechal JD. GaudiMM: A modular multi‐objective platform for molecular modeling. J Comput Chem. 2017;38(24):2118-26. doi: 10.1002/jcc.24847, PMID 28605037.

  82. Sanchez Aparicio JE, Sciortino G, Herrmannsdoerfer DV, Chueca PO, Pedregal JR, Marechal JD. Gpathfinder: identification of ligand-binding pathways by a multi-objective genetic algorithm. Int J Mol Sci. 2019;20(13):3155. doi: 10.3390/ijms20133155, PMID 31261636.

  83. Verdonk ML, Cole JC, Hartshorn MJ, Murray CW, Taylor RD. Improved protein-ligand docking using gold. Proteins. 2003;52(4):609-23. doi: 10.1002/prot.10465, PMID 12910460.

  84. Sanchez Aparicio JE, Sciortino G, Mates Torres E, Lledos A, Marechal JD. Successes and challenges in multiscale modelling of artificial metalloenzymes: the case study of POP-Rh2 cyclopropanase. Faraday Discuss. 2022;234:349-66. doi: 10.1039/d1fd00069a, PMID 35147145.

  85. Siegbahn PE, Blomberg MR. Density functional theory of biologically relevant metal centers. Annu Rev Phys Chem. 1999;50:221-49. doi: 10.1146/annurev.physchem.50.1.221, PMID 15012412.

  86. Bochevarov AD, Harder E, Hughes TF, Greenwood JR, Braden DA, Philipp DM. Jaguar: a high-performance quantum chemistry software program with strengths in life and materials sciences. Int J Quantum Chem. 2013;113(18):2110-42. doi: 10.1002/qua.24481.

  87. Tolbatov I, Re N, Coletti C, Marrone A. An insight on the gold(I) affinity of golB protein via multilevel computational approaches. Inorg Chem. 2019;58(16):11091-9. doi: 10.1021/acs.inorgchem.9b01604, PMID 31353893.

  88. Li YL, Mei Y, Zhang DW, Xie DQ, Zhang JZ. Structure and dynamics of a dizinc metalloprotein: effect of charge transfer and polarization. J Phys Chem B. 2011;115(33):10154-62. doi: 10.1021/jp203505v, PMID 21766867.

  89. Rodzik A, Railean V, Pomastowski P, Zuvela P, Wong MW, Sprynskyy M. Study on silver ions binding to β-lactoglobulin. Biophys Chem. 2022;291:106897. doi: 10.1016/j.bpc.2022.106897, PMID 36240661.

  90. Wang G, Kincaid B, Zhou H, Annaberdiyev A, Bennett MC, Krogel JT. A new generation of effective core potentials from correlated and spin-orbit calculations: selected heavy elements. J Chem Phys. 2022;157(5):054101. doi: 10.1063/5.0087300, PMID 35933201.

  91. Yang ZZ, Wang JJ, Zhao DX. Valence state parameters of all transition metal atoms in metalloproteins development of abeemσπ fluctuating charge force field. J Comput Chem. 2014;35(23):1690-706. doi: 10.1002/jcc.23676, PMID 25042901.

  92. Pomastowski P, Sprynskyy M, Zuvela P, Rafinska K, Milanowski M, Liu JJ. Silver lactoferrin nano complexes as a potent antimicrobial agent. J Am Chem Soc. 2016;138(25):7899-909. doi: 10.1021/jacs.6b02699, PMID 27263865.

  93. Joy A, Biswas R. Role of metal cofactor in enhanced thermal stability of azurin. J Phys Chem B. 2023;127(20):4374-85. doi: 10.1021/acs.jpcb.3c00318, PMID 37183371.

  94. Sheikhzadeh A, Safaei M, Fadaei Naeini V, Baghani M, Foroutan M, Baniassadi M. Multiscale modeling of unfolding and bond dissociation of rubredoxin metalloprotein. J Mol Graph Model. 2024;129:108749. doi: 10.1016/j.jmgm.2024.108749, PMID 38442439.

  95. Bhattarai A, Miao Y. Gaussian accelerated molecular dynamics for elucidation of drug pathways. Expert Opin Drug Discov. 2018;13(11):1055-65. doi: 10.1080/17460441.2018.1538207, PMID 30371112.