Assessment of Molecular Mechanics-based Zn2+ Models in Mono- and Bimetallic Ligand Binding Sites

Zn2+ ions play an important role in biology, but accurate sampling of metalloproteins using Molecular Mechanics remains challenging. Several models have been proposed to describe Zn2+ in biomolecular simulations, ranging from nonbonded models, employing classical 12-6 Lennard-Jones (LJ) potentials or extended LJ-potentials, to dummy-atom models and bonded models. We evaluated the performance of a large variety of these Zn2+ models in two challenging environments for which little is known about the performance of these methods, namely in a monometallic (Carbonic Anhydrase II) and a bimetallic ligand binding site (metallo-β-lactamase VIM-2). We focused on properties which are important for a stable, correct binding site description during molecular dynamics (MD) simulations, because a proper treatment of the metal coordination and forces are here essential. We observed that the strongest difference in performance of these Zn2+ models can be found in the description of interactions between Zn2+ and non-charged ligating atoms, such as the imidazole nitrogen in histidine residues. We further show that the nonbonded (12-6 LJ) models struggle most in the description of Zn2+-biomolecule interactions, while the inclusion of ion-induced dipole effects strongly improves the description between Zn2+ and non-charged ligating atoms. The octahedral dummy-atom models result in highly stable simulations and correct Zn2+ coordination, and are therefore highly suitable for binding sites containing an octahedral coordinated Zn2+ ion. The results from this evaluation study in ligand binding sites can guide structural studies of Zn2+ containing proteins, such as MD-refinement of docked ligand poses and long-term MD simulations.

Introduction modelled explicitly. However, this model generally requires system-specific parameterization of the newly introduced bonds and angles (dihedral angles are generally set to zero in this approach) as well as refitting partial charges of the metal ion(s) and the coordinating residues. Multiple strategies have been proposed to parameterize the bonded parameters. Thereof, the Seminario 19 method is one of the most regularly applied methods, in which the force constants are derived from the Cartesian Hessian matrix. However, this method requires Quantum Mechanics (QM) calculations, which, besides the requirement for access to suitable QM-software and resources and thereby limiting the set-up of high-throughput calculations, may introduce problems. For example, the metal-complex could be unstable in the QM-geometry optimization, which makes the parameterization even more difficult. To avoid these drawbacks, Yu et al. 20 recently presented the Extended Zinc Amber Force Field (EZAFF), which allows a fully empirical derivation of the bonded parameters required in the bonded model.
The major disadvantage of the bonded model is the lack of ligand exchange during the simulation. Since all metal ligands are explicitly bonded to the metal ion, a change of the metal environment cannot be sampled within this approach. Already in 1990, Åqvist and Warshel 21 described a dummy-atom model (also called multisite model) for Mn 2+ , whereas in 1999, Pang 22 was the first to parameterize this type of model for tetrahedral Zn 2+ . In this method, non-interacting dummy atoms carrying a portion of the metal's mass and charge are placed around the metal ion, in a tetrahedral or octahedral geometry. The dummy-atom models can rotate freely around the metal ion without any energy penalty, but bonds between the metal ion and the dummy atoms keep them in the correct distance and geometry around the metal ion. Therefore, this alternative approach allows exchange of metal-ligating atoms, but still enforces/supports a certain coordination geometry due to the charge delocalization in the dummy atoms. More recently, Duarte et al. 23 (re)parameterized the dummy-atom model for octahedral Zn 2+ for the OPLS-AA force field, with the main focus to reproduce metal-oxygen distances, solvation free energy and coordination numbers. Two years later, Jiang et al. 24 further refined the dummy-atom model parameters, using revised reference solvation free energies with Tissandier's 25 proton hydration free energies.
The presence of all those different models illustrates the complexity of a classical description of a metal ion. Since the nonbonded models, including the dummy-atom models, are parameterized in bulk water, there is little information about their performance in a protein environment. Yu et al. 20 compared their EZAFF-derived bonded model in several Zn 2+ -containing systems with their bonded and nonbonded models and a selection of popular semi-empirical QM methods, based on geometry-optimized structures. However, in MD simulations, it is of great value to simulate around an equilibrium state, requiring analysis of the preferred conformational state of these Zn 2+ models. Therefore, we performed long time scale MD simulations of the above mentioned models, analyzing the strengths and weaknesses of these Zn 2+ models to reproduce experimental protein conformations around Zn 2+ ions. This systematic evaluation was performed for Zn 2+ ions in ligand binding sites, since proper modeling of the geometry of the zinc's environment is most challenging in these cases, mostly because of the rather high amount of flexibility and presence of multiple potential Zn 2+ ligating atoms. Additionally, a large amount of metalloproteins, such as metallo-βlactamases, contain two Zn 2+ ions in the binding site. Therefore, the performance of these Zn 2+ models was evaluated in a monometallic as well as a bimetallic ligand binding site, using Carbonic Anhydrase II and the Verona integrin-encoded metallo-β-lactamase (VIM-2) as model systems. Figure 1. Evaluation systems used in this study. Carbonic Anhydrase II (CAII) was applied as model system for monometallic systems, while the β-lactamase VIM-2 was applied as model system for bimetallic systems. The proteins are shown in cartoon representations. The Zn 2+ binding site, containing all coordinating residues and the inhibitor, as well as water molecules and protein residues which are close to the binding site and contain a heavy atom (potentially) able to coordinate the Zn 2+ ion(s), are shown as licorice. The Zn 2+ ions are illustrated as grey spheres. Dotted lines represent coordinate bonds.

Materials and Methods
Force Fields and Potential Function Forms. The nonbonded model in the AMBER force field is described in the following form: where rij is the distance between particle i and j, and Rmin,ij and εij are the distance between this particle pair at which the Lennard Jones (LJ) potential reaches its minimum and the LJ-well depth respectively, applying the Lorentz-Berthelot combining rules. qi and qj are the partial charges of the respective particles.
To improve the description of highly charged systems, Li and Merz 11 presented a modified LJ potential for all interactions with metal ions, including a 1/r 4 term describing the charge-induced dipole interactions. This potential is known as the 12-6-4 LJ-type potential because of its mathematical form: introducing κ as a scaling factor with the unit Å -2 . 4 is defined as 4 = 6 for interactions between the metal ion and the oxygen of water molecules, which is parameterized for a large set of divalent ions by Li and Merz. 11 For the interactions of other atom types with metal ions, 4 is defined as follows: where α0 is the polarizability of the metal ion or respective atom type.
Zn 2+ Models. We evaluated four classical nonbonded (NB) parameter sets: the parameters from Merz 26 , which is the standard parameter set in the AMBER ff03 force field, as well as the three parameter sets from Li et al. 10 The latter were optimized to reproduce the hydration free energy (NB-Li HFE ), ion-oxygen distance (NB-Li IOD ) or a compromise of both (NB-Li CM ), where NB-Li CM is the default parameter set in AMBER ff14SB. The major difference between these parameter sets lies in the value of the LJ well-depth, which ranges from 7.2·10 -4 kcal·mol -1 to 1.5·10 -2 kcal·mol -1 . Because of the recent GPU implementation of the 12-6-4 LJ-type potential in the Amber MD-code, the 12-6-4 LJ-type parameter set could be evaluated in this study as well. Additionally, the three available Zn 2+ dummy atom models (DU) were evaluated, containing one tetrahedral and two octahedral models. The dummy-atoms for the dummy-atom models were placed around the Zn 2+ ions at a distance matching the equilibrium distance of the dummy-atom model applied. For the tetrahedral dummy atom model from Pang 27 , the revised force constant from Park et al. 28 was used for the bonds involving dummy atoms. For the DU-Pang mod model, the binding site was prepared according to the simulation setup from Pang 27-28 : metal-coordinating histidine residues were treated in their deprotonated form (histidine anion), and in the bimetallic binding site of VIM-2, Glu118 was protonated to form a hydrogen bond with the hydroxide ion. Additionally, all residues forming a hydrogen bond with a histidine anion were protonated. The parameters for the histidine anion and hydroxide ion were taken from Pang (2001) 27 .
Last, the performance of the bonded model (BM) was analyzed, with parameters derived from either the Seminario method or EZAFF. Since the partial charges used in this model are both retrieved from the same RESP procedure 29 , these values are identical between these two sets. MCPB.py 17 was used to parameterize the divalent bonds for the bonded models: the positions of the hydrogen atoms from the large model generated by MCPB.py were optimized at the B3LYP/6-31G(d) level of theory. This optimized geometry was used for the RESP fitting, applying the ChgModB scheme, i.e. restraining all charges of the backbone heavy atoms to the values from the ff14SB force field, since this gave the best results in a study by Peters et al 18 . The van der Waals radius of the Zn 2+ ion was set to 1.395 Å, taken from the IOD set of Li et al 10 . For the empirical model (BM-EZAFF), the bonded parameters were derived from the Extended Zinc Amber Force Field (EZAFF) 20 . For the Seminario model (BM-Seminario), a small model containing the Zn 2+ ion and the coordination residues and ligand was generated with MCPB.py. The geometry of this small model was optimized at the B3LYP/6-31G(d) level of theory. The bonded parameters were derived applying the Seminario method on this geometry optimized small model. For VIM-2, the small model was initially optimized at the TPSS/def2-TZVP level of theory prior to the application of the Seminario method, to improve the description of the complex electronic structure in the bimetallic binding site. Additionally, a small van der Waals radius was placed on the hydrogen of the hydroxide ion, taken from the GAFF2 force field. All dihedrals involving a Zn 2+ ion were set to zero, in agreement with the strategy applied in the MCPB procedure. The most relevant parameters of the evaluated Zn 2+ models are shown in Table 1.
Simulation Setup. Carbonic Anhydrase II (CAII) in complex with the RA1 inhibitor (PDB ID: 5NXG) resolved at 1.2 Å resolution was used for the evaluation of Zn 2+ models in monometallic systems, while the metallo-β-lactamase VIM-2 with the ANT-431 inhibitor (PDB ID: 6HF5), resolved at 1.8 Å resolution was applied as model system for the evaluation in bimetallic binding sites. [30][31] During the ligand parameterization procedure, the ligands were protonated with Open Babel 32 (v. 2.3.2.) at pH 7.0, keeping a deprotonated NHfrom the sulfonamide group of the CAII inhibitor, as this has been shown to be the state the ligand binds to the Zn 2+ ion. 11,[32][33] For the bonded and van der Waals parameters of the ligand atoms, the General Amber Force Field 34 (GAFF) parameters were used. Charges of the ligand atoms for all but the bonded models were derived applying the RESP procedure based on a Merz-Singh-Kollman population analysis 35-36 at a QM-geometry optimized structure of the inhibitor, which was performed at the HF/6-31G(d) level of theory, with Gaussian09 (revision E.01) 37 .
The protein-ligand complexes were prepared as follows: from all multi-resolved residues, the first occurrence was selected, and all water molecules within a sphere of 8 Å of the Zn 2+ ion(s) were preserved. The protonation states of ionizable protein residues were predicted with the PROPKA3.0 software package 38-39 at pH 7.0, followed by a visual check. For all standard protein residues, the AMBER ff14SB 40 force field was applied. The system was solvated in a rectangular box containing TIP3P 41 water, applying a buffer region of 12 Å around the protein atoms. The system was neutralized with Na + ions using parameters from Joung and Cheatham 42 , using the LEaP module of the Amber18/AmberTools18 software package 16 . The bridging hydroxide ion in VIM-2 was parameterized as described in a previous study by Marion et al. 43 containing a -1.3 e and +0.3 e partial charge for the oxygen and hydrogen atom respectively. The polarizability of the hydroxide oxygen, needed to calculate the C4 parameter in the 12-6-4 LJ-type model, was set to 2.03 Å 3 , taken from Shannon et al. 44 , while the polarizability of the hydrogen atom was set to zero, analogously to the hydrogen of water. The Zn 2+ parameters from Table 1 were applied, and for the 12-6-4 LJ-type model, the C4 parameters were added with ParmEd.
Simulation Procedure. An energy minimization was performed using sander from Amber18/AmberTools18 with the XMIN method (ntmin=3), thereby gradually adjusting the box size to bring the density from 0.8 g·cm -3 to 1.0 g·cm -3 , in steps of 0.02 g·cm -3 . During this minimization, a 20.0 kcal·mol·Å -2 positional restraint was applied to all protein atoms. When the target density was reached, an additional minimization was performed, applying the positional restraint solely to the binding site, defined as all residues with at least one atom in a sphere of 5 Å around the Zn 2+ ion(s). During the heat-up procedure, the temperature was gradually increased to 300 K, while decreasing the positional restraint. The precise methodology is provided in the Supporting Information (Table S2 in the Supporting Information). The Langevin thermostat 45 was applied with a collision frequency of 4.0 ps -1 . During the MD simulations in the NPT ensemble, the Berendsen barostat 46 was applied to keep the pressure at 1 bar, with a relaxation time of 1 ps and compressibility of 44.6·10 -6 ·bar -1 . Periodic boundary conditions were applied and the SHAKE algorithm 47 was used to constrain all bonds involving a hydrogen atom at their equilibrium distance. A cut-off distance of 12 Å was used for all nonbonded interactions, while the particle mesh Ewald method 48 was applied to describe long range electrostatics. The simulations were performed with an integration step size of 1 fs. After the heat-up, a 200 ns production MD simulation was performed at 300 K in a NPT ensemble, saving the atom coordinates and velocities every 10 ps. The heat-up and production simulations were conducted in three replicas, with the pmemd.cuda MD engine from Amber18.
Trajectory analyses. From each replicon, the first 100 ns were considered as equilibration, thus the last 100 ns from the three replicas were merged, and analyzed as a single trajectory of 30000 frames. The distances and distance RMSDs between Zn 2+ and the ligating atoms were calculated with cpptraj 49 . The coordination geometry was determined every 2 ns for each Zn 2+ individually, using a python package FindGeo 50 . Every non-carbon heavy atom within 2.8 Å of the Zn 2+ ion was considered as ligating atom, except the sulfur from the sulfonamide moiety in the CAII ligand. Rmin/2, ε and κ are provided in Å, kcal·mol -1 , and Å -2 respectively. a Dummy atom models (DU) additionally differ in Zn 2+ -dummy atom bond/angle/dihedral parameters. b Bonded models (BM) differ in parameters regarding Zn 2+ -bonded residues, provided in the Supporting Information (Table S1). c Retrieved from RESP fit.

Results
We evaluated a large set of available Zn 2+ models, which can be subdivided in four model types: classical nonbonded models (NB), the 12-6-4 LJ-type model (LJ1264), dummy-atom models (DU) and bonded models (BM). The performance of these methods was evaluated in the binding sites of CAII and VIM-2 ( Figure 1), with the protonation state of the ligands and the binding site residues as supported by experiment (see Methods section). For the dummy atom model from Pang 27-28 , we conducted the performance evaluation with two different protonation states of the coordinating residues. The first evaluated setup (DU-Pang) is the protonation state identical as used for the other models while the second setup (DU-Pang mod ) is simulated with deprotonated histidine residues and a protonated Asp118. The latter is the simulation setup equivalent to the setup used by Pang. Figure 2. RMSD of interatomic distances (dRMSD) between Zn 2+ and selected binding site atoms with respect to the average simulated interatomic distance of the respective atom pair, as a measure for the stability of the sampled Zn 2+ binding site. The distances to histidine residues are defined as the distance between the respective Zn 2+ ion to the ligating nitrogen of the histidine. The distances to an aspartate or glutamate residue are measured between the respective Zn 2+ ion and the carboxyl carbon, since these side chains can rotate freely and contain two chemically equivalent oxygen atoms. A: distance RMSD in CAII. O1 (ligand) and O2 (ligand) represent the sulfonamide oxygen atoms, which are chemically equivalent, but an exchange of their position will lead to large changes in the binding pose of the ligand. N (ligand) represents the amine nitrogen of the sulfonamide moiety of the ligand. B: distance RMSD in VIM-2. The Zn2 2+ -CO2 (ligand) distance represents the distance between Zn2 2+ and the carboxyl carbon from the inhibitor, where Zn2 2+ -N (ligand) is measured towards the nitrogen of the ligand's thiazole moiety To analyze the integrity of the Zn 2+ binding site, the Root-Mean-Square Deviation was calculated over the interatomic distances (dRMSD) between the Zn 2+ ion and a selection of atoms defining the structure of the binding site ( Figure 2). The dRMSD is calculated with respect to the average distance (over all three replicas) simulated between the respective atom pair. Small dRMSD values thus indicates a stable ligation of the Zn 2+ ion, because the interatomic distance between the Zn 2+ ion and a ligating atom only slightly deviates from the average distance observed in the respective simulations. Large dRMSD values on the contrary indicate unstable ligation of the Zn 2+ ion. The simulated Zn 2+ binding sites by the nonbonded models show the largest distortion of the binding sites of all evaluated Zn 2+ models, with the most unstable MD trajectories for NB-Merz and NB-Li CM . Besides the high average dRMSD values in the nonbonded models, the large error bars, especially for measures including histidine residues, illustrate that no stable conformation could be sampled within all replicate simulations. Addition of the ion-induced dipole effects via the 12-6-4 LJ-type potential improves the reproduction of the binding site significantly, especially for the Zn 2+ -His coordination bonds. Regarding the DU-Pang models, only simulations applying DU-Pang mod show stable tetrahedral Zn 2+ ions (CAII, and Zn1 2+ in VIM-2). Remarkable is the strong performance of DU-Duarte and DU-Jiang, as these dummy-atom models simulate the VIM-2 binding site even closer to the experimental structure as the bonded models, illustrated by the dRMSD < 0.5 Å for all Zn 2+ -ligating atom distances. These dummy-atom models did not only reproduce the binding site for the octahedral coordinated Zn 2+ ions, but also for the binding site around tetrahedral Zn 2+ .  a Native coordination geometry, as observed in the X-Ray structure for CAII (PDB ID: 5NXG). b All coordination geometries which were visited less than 5 % in all evaluated Zn 2+ models in both CAII and VIM-2 were merged in this category. A complete list of all sampled coordination geometries for CAII is provided in the Supporting Information (Table S3).
Another important factor is the coordination geometry sampled by the Zn 2+ models. This coordination geometry determines the exact position of the Zn 2+ ligating atoms, which plays an important role in either stabilization of the protein structure, or introducing strain on certain bonds, for example to catalyze a chemical reaction. Therefore, we followed the coordination geometry of the Zn 2+ ion(s) sampled during the simulations by determining the coordination geometry every 2 ns simulation time. The coordination geometry of the Zn 2+ ion observed in the X-Ray structure CAII is tetrahedral, but this is only the most sampled coordination geometry by NB-Li HFE , DU-Pang mod , and the bonded models (Table 2). Almost all remaining Zn 2+ models prefer an octahedral Zn 2+ geometry, with 83% to 96% occurrence. The trigonal bipyramidal and square pyramidal geometries were barely sampled.   (Table S4). c All coordinating ligands were water molecules or a hydroxide ion, and in one replicon the carboxylate oxygen of the ligand.
A similar trend can be observed in the bimetallic system, VIM-2 ( Table 3). The octahedral coordination geometry for Zn1 2+ is clearly preferred by most nonbonded models, while a tetrahedral coordination geometry is found in the X-Ray structure of VIM-2. Besides the bonded models, only with the NB-Li HFE and DU-Pang mod models a tetrahedral coordination is sampled predominantly for Zn1 2+ . Strikingly, the preferred coordination geometry for Zn2 2+ , which is octahedral according to the X-Ray structure, was not always sampled in this octahedral geometry. With the NB-Merz, NB-Li HFE , NB-Li CM , and DU-Pang mod models Zn2 2+ is sampled mostly in a tetrahedral geometry. NB-Li IOD and the 12-6-4 LJ-type model, however, simulated Zn2 2+ in an octahedral geometry. We also analyzed the preference of the Zn 2+ models for certain ligating residues during the simulations (Figure 3). We defined a residue as ligating whenever their ligating atom is present within the coordination sphere of the Zn 2+ ion, i.e. a sphere with a 2.8 Å radius around the Zn 2+ ion. A normalized/average coordination number (CN) was calculated, where CN = 1 indicates a stable monodentate coordination between this residue and the Zn 2+ ion. CN < 1 indicates a less stable interaction, since the respective coordination interaction was not constantly present over the entire simulation. CN values > 1 for a protein residue indicates that on average more than one ligating atom contributes to the Zn 2+ coordination, which could occur for bidentate residues having more than one atom able to ligate the Zn 2+ ion (e.g. Asp and Glu). CN > 1 for water molecules indicate that multiple water molecules coordinate the Zn 2+ ion. This analysis shows a clear difference in the preference of the evaluated Zn 2+ models towards specific ligating atoms. For example, the interaction with histidine residues is often not properly sampled, because the coordination is taken over by a charged residue or a water molecule. This evaluation gave further insights in the preference of some Zn 2+ models for charged ligating atoms, or hard or soft bases, which will be discussed below.

Discussion
The performance of the Zn 2+ models was analyzed in the binding sites of Carbonic anhydrase II and β-lactamase VIM-2, as model systems for monometallic and bimetallic binding sites respectively (Figure 1). CAII is very well suitable for the analysis of monometallic Zn 2+ binding sites, since the Zn 2+ ion in CAII is coordinated in a tetrahedral geometry, which is the dominant coordination geometry of Zn 2+ monometallic binding sites. 8 The Zn 2+ ion is coordinated by three histidine residues and the nitrogen of the sulfonamide moiety of the inhibitor, again typical ligating residues for Zn 2+ . Due to the opposite charge between Glu106 and the Zn 2+ ion, we hypothesized that certain Zn 2+ models may overestimate this electrostatic interaction, and thereby falsely consider Glu106 as a coordinating residue. Therefore, the non-coordinating glutamate residue Glu106 close to the Zn 2+ ion (4.0 Å and 5.5 Å between Zn 2+ and both glutamate oxygen atoms respectively) was included in the analysis as well. Furthermore, Glu106 plays an important role in an elaborate hydrogen-bonding network in CAII: this glutamate orients binding site residues such that they can interact with the ligand, an effect which has been described in a variety of studies, illustrating the importance of proper sampling of this residue. [51][52][53] Glu106 interacts with hydrogen bonds via a water molecule to Tyr7, and directly to the side chain of Thr199, which in turn interacts both via the backbone nitrogen and the side chain hydroxyl group to the sulfonamide moiety of the inhibitors.
The binding site of VIM-2 contains two Zn 2+ ions, bridged by a hydroxide ion. 54-56 Zn1 2+ is coordinated in a tetrahedral geometry by three histidine residues and the bridging hydroxide, while Zn2 2+ adopts an octahedral coordination geometry via an aspartate, cysteine and histidine residue, as well as the bridging hydroxide, the thiazole nitrogen and carboxylate oxygen of the inhibitor. Besides the interaction of the ligand with the Zn 2+ ion, the ligand forms several hydrogen bonds with the protein: Arg228 forms a salt-bridge with the carboxylate group of the ligand, and Tyr67 forms an aromatic interaction with the ligand's pyridine moiety. 31 The two differently coordinated Zn 2+ ions make this an interesting model system, allowing evaluation of the Zn 2+ models to discriminate between different coordination geometries. Furthermore, the effect of the presence of a second Zn 2+ ion within a single system on the Zn 2+ model's performance can be studied with this system.
In this study, we evaluated the performance of a variety of Zn 2+ models in these two ligand binding sites, more precisely, in their ability to describe the structure of the Zn 2+ coordination environment properly during MD simulations. We performed this study in ligand binding sites, since we think this is one of the most challenging case for the Zn 2+ models, mainly due to the commonly rather high flexibility of ligand binding sites. Furthermore, ligand binding sites often contain a high number of solvent molecules, which are all able to coordinate Zn 2+ . Another reason to evaluate the performance in ligand binding sites is that proper structural sampling is highly important in the ligand binding site for a variety of practical applications, like (MD-based) molecular docking simulations and refinements, or studies with special focus on resolving chemical reaction pathways.
Classical nonbonded models. Proper sampling of the metal's coordination geometry is very important, mainly because the ion's ability to bind its environment in a certain coordination geometry is one of the most important properties of transition metal ions in proteins. However, this property is often neglected. [57][58] We found that the coordination geometry is least accurately reproduced by the nonbonded models compared to the other models in this study (Table 2 and 3), and that often incorrect ligating atoms are found (Figure 3 and 5), implicating that the standard 12-6 LJ potential is insufficient for proper description of Zn 2+ ions. This observation follows previous findings from a recent evaluation study of Mg 2+ models by Zuo and Liu. 58 At first sight, the NB-Li HFE model seems to work well, since a tetrahedral coordination geometry was sampled for the tetrahedral Zn 2+ ions (CAII, and Zn1 2+ in VIM-2). However, this tetrahedral coordination was not sampled with the correct ligating residues: in all replicas, His119 was replaced by a ligating water molecule and the ligating role of His94 was taken over by Glu106 in one of the replicate simulations ( Figure 3). This artificial solvation of the Zn 2+ ion in CAII was also observed for the NB-Merz and NB-Li CM models, for which most coordination bonds with the histidine residues were broken and replaced by water molecules. Simulations of VIM-2 applying these nonbonded models again showed that all coordination bonds between the Zn 2+ ions and the histidine residues, as well as the ligand's thiazole nitrogen were broken and replaced by oxygen atoms, mostly from water molecules ( Figure 4). Our data reveals a preference of the nonbonded models, with the exception of NB-Li IOD , towards hard bases as ligating atoms (as in water). The only exception is the interaction between Zn2 2+ and Cys198 in VIM-2. This interaction remained intact in all Zn 2+ models, despite the rather soft base character of cysteine. However, the interaction between Cys198 and Zn2 2+ may be prevailed by the coulomb term here, because this residue was modelled in the deprotonated state, which may explain the observed stable interaction. We think that this preference of hard bases over soft bases as ligating atoms by the nonbonded models may partly be a result of the parameterization procedure, which is traditionally performed for a single ion in bulk solvent (water).
Both NB-Merz and NB-Li CM models show a highly similar pattern in dRMSD measures and preferred coordinating residues (Figure 2, 5 and 6). Furthermore, both models preferred an octahedral coordination for CAII and Zn1 2+ in VIM-2, while they simulated a tetrahedral coordination geometry for Zn2 2+ (Table 2 and 3). The instable simulation trajectories and incorrect Zn 2+ coordination sampled by these models result mainly from unstable binding of the Zn 2+ ions, as the Zn 2+ ions move up to 3 Å away from their crystallized position during the simulations ( Figure  S1 and S2 in the Supporting Information). NB-Li IOD , which was specifically optimized to reproduce ion-oxygen distances, properly described the Zn 2+ -His interactions for His96 and His119 in CAII and His116 and His179 in VIM-2, in contrast to the other 12-6 LJ nonbonded models as described above. However, this nonbonded model still failed to predict the interaction with His94 in CAII (in one replicon) and His114 in VIM-2 (in all replicas), as well as the thiazole nitrogen of the ligand in VIM-2, where the latter results in large deviations in the ligand binding pose ( Figure  3 and 4). Additionally, an octahedral coordination geometry was sampled with the NB-Li IOD model for all evaluated Zn 2+ ions (Table 2 and 3), thus this model seems not to be suitable to reproduce tetrahedral binding sites. As described before, this preference for an octahedral coordination geometry could possibly be an artifact from the parameterization in water, since Zn 2+ coordinates water in an octahedral geometry. [59][60]  Influence of ion-induced dipole effects. The idea from Li and Merz to include ion-induced dipole effects in the LJ-equation for all interactions containing a metal ion, i.e. the 12-6-4 LJ-type potential, leads to much better description of interactions between Zn 2+ and particularly soft/borderline bases (Figure 3 and 5). However, we observed that the 12-6-4 LJ-type model almost solely models the Zn 2+ ion in an octahedral coordination geometry (Table 2 and 3). This can be problematic whenever an accurate reproduction of the metal coordination is required, for example if MD simulations are used to refine docked ligand poses, because Zn 2+ often has a tetrahedral coordination in protein structures. 8 To accomplish this octahedral coordination geometry in the tetrahedral binding sites CAII, and Zn1 2+ in VIM-2, additional water molecules were ligated in simulations applying the 12-6-4 LJ-type potential (Figure 3), an effect which was also observed in simulations with other metal ions. 61 Furthermore, our simulations of CAII applying the 12-6-4 LJtype model still overshoot the interaction between Zn 2+ and the charged Glu106, as illustrated by the too short distance between them ( Figure 5). Despite these inaccuracies, this modified LJpotential leads to much better description of Zn 2+ coordination bonds in a protein environment, and thereby increases the stability of the binding site geometries considerably compared to the classical 12-6 LJ nonbonded models (Figure 2). This better metal description also improved the sampling of the ligand binding pose in VIM-2 ( Figure 4): the heavy-atom RMSD of the ligand is 1.8 Å in the simulations applying the 12-6-4 LJ-type potential, using the crystal structure as reference, while the ligand RMSD of the other nonbonded models are higher than 2.8 Å.
Dummy atom models and the influence of protonation states. The simulation of the tetrahedral binding sites can be improved substantially by the DU-Pang mod dummy atom model. This model sampled the tetrahedral binding sites almost exclusively in their correct coordination geometry: 96.7% and 99.3% for CAII, and Zn1 2+ in VIM-2 respectively (Table 2 and 3). However, this result could only be reached applying the (artificial) deprotonated histidine residues and modified protonation states of the other coordinating residues, as used by Pang in his publication as well. 22 The performance worsened drastically once the same protonation states were used as in all other Zn 2+ models, i.e. neutral histidine residues, (named DU-Pang in this study): no stable coordination geometries could be found in these cases, only 20.7% and 4.7% of the simulation showed a tetrahedral coordination for CAII and Zn1 2+ in VIM-2, respectively. Furthermore, following Pang's rules to determine which residues should be protonated, the required protonation state of Asp118 in VIM-2 remains ambiguous. Asp118 should be deprotonated according to Pang's rules because of its interaction with Zn2 2+ , yet it must be protonated because of its hydrogenbonding with the hydroxide ion. Therefore, we performed the simulations applying the DU-Pang mod model with both protonated and deprotonated Asp118: the simulation containing a deprotonated Asp118 resulted in a highly distorted binding site, thereby losing ligand binding (data not shown), while the simulation with a protonated Asp118 stayed rather close to the experimentally observed structure. Therefore, the simulation with the protonated Asp118 represents the DU-Pang mod model in this study. This protonated Asp118 did however not form a hydrogen bond with the hydroxide ion during the simulation, which was the rationalization to protonate this residue. In contrast, Asp118 formed a hydrogen bond with solvent molecules and thereby losing its coordination of the Zn 2+ ion (Figure 3 and 5). Both DU-Pang and DU-Pang mod were not able to reproduce the octahedral binding geometry from Zn2 2+ in VIM-2. With the DU-Pang model an octahedral geometry was simulated in only 0.7% of the simulation time, and with the DU-Pang mod model a tetrahedral geometry was found instead (Table 3). Thus, DU-Pang mod performs very well for tetrahedral Zn 2+ ions, but the performance is highly sensitive on the choice of protonation state of the coordinating residues and is not suitable to reproduce another coordination geometry than tetrahedral. Figure 5. Interatomic distance between Zn 2+ and selected binding site atoms, as a measure for binding site integrity for the Zn 2+ ion in CAII (A), and Zn1 2+ (B) and Zn2 2+ (C) in VIM-2. The bars represent the average distance between the Zn 2+ ion and the ligating residue during the simulation, while the dotted lines represent the value of the respective distance in the X-Ray structure. All distances to an aspartate or glutamate residue are measured to the side-chain's carboxyl carbon, since both oxygen atoms are equivalent and can freely take over each-other's function.
The other Zn 2+ dummy atom models, DU-Duarte and DU-Jiang, the latter being a reparameterization of DU-Duarte based on revised solvation free energies 24 , were both parameterized for octahedral Zn 2+ ions, and therefore contain six dummy-atoms. [23][24] Despite their coordination geometry specific parameterization, they also performed quite well for the tetrahedral Zn 2+ ions in this study. The six dummy-atoms enforce an octahedral coordination geometry, also for the tetrahedral Zn 2+ ions (Table 2 and 3), but this had only a minor effect on the positions and stability of the coordinating residues ( Figure 5). The additional coordination partners required to fulfill the octahedral coordination geometry for Zn1 2+ in VIM-2 were mainly water molecules, which we also observed for simulations applying the 12-6-4 LJ-type potential (Figure 3). In CAII, in both, the DU-Duarte and DU-Jiang models, a sulfonamide oxygen was used as an additional ligating atom, but this barely affected the ligand binding orientation. Using the DU-Duarte model, an additional ligation with Glu106 was observed in all replicas, yet this glutamate kept its important interaction with Thr199, which was shown to be a conserved interaction in CAII. With the DU-Jiang model, ligation of Glu106 was found only in one of the replicas, while a water molecule was used in the other two. Thus, this resulted in an almost perfect reproduction of the Zn 2+ binding site and in very stable simulation trajectories (Figure 2), both for the monometallic and bimetallic system, despite the wrong coordination environment of the tetrahedral Zn 2+ ions. Both models, DU-Duarte and DU-Jiang, perform very well for the octahedral Zn2 2+ in VIM-2, applying the correct ligating residues over the entire simulation ( Figure 3).
Finally, we also analyzed if a tetrahedral and octahedral dummy-atom model can be combined in the bimetallic system VIM-2. Therefore, we performed simulations modelling the tetrahedral Zn1 2+ by either DU-Pang or DU-Pang mod , and the octahedral Zn2 2+ by DU-Jiang. The simulations applying the combined DU-Pang mod and DU-Jiang models resulted in highly stable simulations and reproduced the coordination geometry of both ions correctly ( Figure S3 and Table S4 in Supporting Information). However, in the combined model simulations, one still needs to (artificially) deprotonate the residues around Zn1 2+ in order to sample a correct coordination geometry. Simulations with protonated and deprotonated Asp118 were similar, whereas the first performed best in the reproduction of the metal environment.
Restraining metal-coordination via explicit bonds. The bonded models reproduced the Zn 2+ sites rather accurately, which could be expected because of the system-specific parameterization and application of explicit bonds between the Zn 2+ ions and ligating atoms. Correct usage of the ligating residues, proper reproduction of the coordination geometry and highly stable simulation trajectories were observed. It is however surprising that simulations of VIM-2 applying the dummy-atom models DU-Duarte and DU-Jiang were even more stable than the bonded models, and the simulated bimetallic binding site geometries were closer to the experimental structure than for both bonded models, while not having the disadvantage of fixed ligating residues. Especially the distance between the two Zn 2+ ions (which were connected via explicit bonds in the bonded model via the hydroxide ion) was simulated closer to the experimental value by these dummy-atom models, as shown in Figure 5. It needs to be noted here that BM-Seminario is parameterized based on a QM-geometry optimized structure, during which the distance between the Zn 2+ ions reduced to 3.7 Å, with respect to 4.0 Å in the X-Ray structure. However, the distance between the two Zn 2+ ions during the BM-Seminario simulations was even smaller (3.3 Å). Finally, we barely observed any performance differences between BM-Seminario and BM-EZAFF, indicating that the empirical derivation of the bonded parameters works well, both for monometallic and bimetallic binding sites. Figure 6. Interatomic distance between Zn 2+ and electronegative atoms from the ligand in CAII (A) and VIM-2 (B). The bars represent the average distance between the Zn 2+ ion and the ligating atom during the simulation, while the dotted lines represent the value of the respective distance in the X-Ray structure. Zn2 2+ -CO2 in VIM-2 is measured from the Zn 2+ ion to the carboxyl carbon of the ligand, since both oxygen atoms are equivalent and can freely take over each-other's function.

Conclusions
Based on our evaluation of a large variety of Zn 2+ models, we show that the nonbonded models applying the standard 12-6 LJ potential are too limited to accurately describe the conformation of a flexible Zn 2+ -containing protein binding site. Especially interactions between the Zn 2+ ion and non-charged ligating atoms are underestimated, such as the imidazole nitrogen in histidine residues. These ligating atoms are often replaced by water molecules or charged residues during the simulation. The performance can be improved by the introduction of the ion-induced dipole interactions via the 12-6-4 LJ-type potential, which thanks to the GPU-implementation in AMBER18 is now also applicable for long time scale simulations. However, we showed that the 12-6-4 LJ-type model strongly prefers an octahedral coordination geometry and overestimates interactions with charged atoms, but this had only a minor effect on the binding site conformations sampled in this study. The dummy-atom models are a good alternative to the nonbonded models, from which DU-Pang mod reproduces a tetrahedral coordination geometry fairly well, but its performance depends strongly on the not always rational protonation state of the metal's environment. Since many Zn 2+ ions in proteins are found in a tetrahedral coordination geometry, the development of a tetrahedral dummy-atom model parameterized specifically in a protein environment would be welcome. The octahedral dummy-atom models DU-Duarte and DU-Jiang are good choices for Zn 2+ ions coordinated in an octahedral geometry because of their solid performance in this study regarding the reproduction of the metal's environment. We also observed that these models perform rather well for tetrahedral coordinated Zn 2+ , but still attract water molecules to fill the coordination sphere. The difference in performance between DU-Duarte and DU-Jiang is minor, whereas the latter showed the best performance in this study. These octahedral dummy-atom models perform very similar, in the case of the bimetallic VIM-2 binding site even slightly better than the bonded models in this study. However, the latter may presumably still be the best choice whenever a fixed coordination geometry needs to be sampled, especially if the required geometry is disturbed and does not match any of the parameterized geometries. The combination of different dummy models based on their coordination geometry in the bimetallic site led to very promising results, indicating that it would be very useful to provide parameter sets for these models for different coordination environments. Because of the almost identical performance of BM-Seminario and BM-EZAFF in this evaluation study, the EZAFF method to parameterize the bonded model will often suffice, having the advantage of being independent of QM-calculations for its parameterization, except for charge fitting.

Supporting Information
The Supporting Information is available free of charge (PDF). Figure S1. Zinc displacement. Figure S2. Sampled coordination geometries in CAII. Figure S3. Performance of combining DU-Pang and DU-Jiang for the bimetallic system VIM-2. Table S1. Bonded parameters applied in bonded models. Table S2. Overview of heat-up protocol applied in this study. Table S3. Full list of sampled Zn 2+ coordination geometries in CAII. Table S4. Full list of sampled Zn 2+ coordination geometries in VIM-2.