Heterogeneous behavior of metalloproteins toward metal ion binding and selectivity: insights from molecular dynamics studies

About one-third of the existing proteins require metal ions as cofactors for their catalytic activities and structural complexities. While many of them bind only to a specific metal, others bind to multiple (different) metal ions. However, the exact mechanism of their metal preference has not been deduced to clarity. In this study, we used molecular dynamics (MD) simulations to investigate whether a cognate metal (bound to the structure) can be replaced with other similar metal ions. We have chosen seven different proteins (phospholipase A2, sucrose phosphatase, pyrazinamidase, cysteine dioxygenase (CDO), plastocyanin, monoclonal anti-CD4 antibody Q425, and synaptotagmin 1 C2B domain) bound to seven different divalent metal ions (Ca2+, Mg2+, Zn2+, Fe2+, Cu2+, Ba2+, and Sr2+, respectively). In total, 49 MD simulations each of 50 ns were performed and each trajectory was analyzed independently. Results demonstrate that in some cases, cognate metal ions can be exchanged with similar metal ions. On the contrary, some proteins show binding affinity specifically to their cognate metal ions. Surprisingly, two proteins CDO and plastocyanin which are known to bind Fe2+ and Cu2+, respectively, do not exhibit binding affinity to any metal ion. Furthermore, the study reveals that in some cases, the active site topology remains rigid even without cognate metals, whereas, some require them for their active site stability. Thus, it will be interesting to experimentally verify the accuracy of these observations obtained computationally. Moreover, the study can help in designing novel active sites for proteins to sequester metal ions particularly of toxic nature.


Introduction
Proteins utilize a large array of cofactors to achieve a variety of structures and functions, among which metal ions are known to perform critical functions and differ from organic cofactors. Approximately, one-third of all proteins are estimated to be metalloproteins (Barondeau & Getzoff, 2004;Waldron & Robinson, 2009). An analysis of the human genome suggests that~3000 (~10% of all encoded) proteins participate in Zn 2+ binding (Andreini, Banci, Bertini, & Rosato, 2006). Apart from playing key role in the metabolism and survival of every organism, metal ions are also involved in microbial virulence (Anderson et al., 2009;Chowdhury, Sahu, & Das, 1996;Diaz-Ochoa, Jellbauer, Klaus, & Raffatellu, 2014;Kim, Watanabe, Shirahata, & Watarai, 2004;Counago, McDevitt, Ween, & Kobe, 2012). The functionality of metalloprotein depends on subtle interactions between the electronic properties of the metal ion, dictated by its coordination chemistry and the stability of protein conformations. Sensor and regulatory metalloproteins are known to energetically match metal ion coordination chemistry with protein conformational change at a binding affinity appropriate for their biological function (Benoit & Maier, 2011;Counago et al., 2014;McDevitt, Ogunniyi, Valkov, Lawrence, & Kobe, 2011). Although metal ions are necessary for cellular functions, they require strict regulation to hinder toxic side reactions or misincorporation into non-native binding sites (Hollenstein et al., 2009;Wang, Lu, Zhou, Wang, & Liu, 2013). Thus, the concentrations of 'free' metal ions are under the control of transcription factors and appear to be exceedingly low in cells (Finney & O'Halloran, 2003). However, the mechanism of protein chain conformational adaptation to the requirements of the metal coordination has not yet been completely understood. The composition of amino acid residues and their separation in the primary structure is insufficient to predict the conformation for metal binding (Harding, 2004). Here, in this study, we have performed molecular dynamics (MD) simulations of seven proteins to investigate whether the cognate metal ion bound to its respective protein is replaceable with other similar metal ions. The proteins (metals) studied in the present work are phospholipase A 2 (Ca 2+ ), sucrose phosphatase (Mg 2+ ), pyrazinamidase (Zn 2+ ), cysteine dioxygenase (Fe 2+ ), plastocyanin (Cu 2+ ), anti-CD4 antibody Q425 (Ba 2+ ), and synaptotagmin 1 C 2 B domain (Sr 2+ ).
The enzyme phospholipase A 2 (PLA 2 ) catalyzes the hydrolysis of the sn-2 ester bond of phospholipids and plays an important role in a number of physiologically significant processes such as inflammation, blood platelet aggregation, and acute hypersensitivity (Pieterson, Vidal, Volwerk, & De Haas, 1974;Thunnissen et al., 1990;White, Scott, Otwinowski, Gelb, & Sigler, 1990). The enzyme has a digestive function and thus is found in large amounts in mammalian pancreatic juice and in snake and bee venoms (Six & Dennis, 2000). A calcium ion (Ca 2+ ) is required for the enzymatic activity and is stabilized by five protein coordination including two carboxyl atoms of Asp49, carbonyl main chain oxygen atoms of Tyr28, Gly30, and Gly32, and two water molecules (Dijkstra, Kalk, Hol, & Drenth, 1981).
The enzyme sucrose phosphatase (SPP) catalyzes the final step in the sucrose biosynthesis by converting sucrose-6-phosphate to sucrose (Hawker & Hatch, 1966;Leloir & Cardini, 1955). In most plants, sucrose is used to store starch, oils and also plays role under drought, salinity, or cold stress. Thus, the enzyme SPP links photosynthetic and non-photosynthetic processes in addition to its adaptive responses to abiotic stresses. The enzyme SPP belongs to the HAD superfamily which are characterized by three conserved motifs ( (Aravind, Galperin, & Koonin, 1998;Collet, Stroobant, Pirard, Delpierre, & Van Schaftingen, 1998;Wang, Kim, Jancarik, Yokota, & Kim, 2001). The residues, in particular the first aspartate, from motif I are essential for the catalytic activity of the enzymes of the HAD superfamily. Many proteins from this superfamily require a divalent metal ion Mg 2+ which is stabilized by residues from motif I and motif III (Morais et al., 2000;Wang et al., 2001). The residues serine or threonine of motif II hydrogen bonds with a phosphoryl oxygen in the substrate anchoring its correct orientation for the nucleophilic attack by the aspartate residue of the motif I .
Cysteine dioxygenase (CDO), a non-heme ironbinding protein, catalyzes the conversion of cysteine to cysteine sulfinate in its ferrous form and plays a crucial role in intracellular sulfur homeostasis (Dominy et al., 2007;Stipanuk, Simmons, Karplus, & Dominy, 2011). Mutations in CDO have been linked with rheumatoid arthritis and neurodegenerative disorders (Stipanuk, 2004). The requirement of iron for the enzymatic activity of CDO has been a matter of discussion. Several studies, in the past, have shown activity of the enzyme without added iron in the pH range of~6 to~9 (Chai, Bruyere, & Maroney, 2006;Gardner, Pierce, Fox, & Brunold, 2010;Karplus & Diederichs, 2012;McCoy et al., 2006;Pierce, Gardner, Bailey, Brunold, & Fox, 2007). However, a recent study suggests that iron in its ferrous form is required for the enzymatic activity in the pH range of~5.5 to~6.5 (Driggers et al., 2013).
Plastocyanin (Pc) is a Cu 2+ -binding protein involved in electron-transfer shuttle from cytochrome f of the cytochrome b 6 f complex of photosystem II to P700 + of photosystem I in chloroplasts of algae, higher plants, and many cyanobacteria (Freeman & Guss, 2001). Plastocyanin bound to Cu 2+ (Cu 2+ Pc) is reduced by cytochrome f to Cu + Pc which diffuses through the lumen until it binds to P700 + where it is oxidized to Cu 2+ Pc releasing an electron (e − ). Algae and cyanobacteria lacking plastocyanin use cytochrome c 6 for this task (Sandmann, Reck, Kessler, & Boger, 1983). Both the proteins plastocyanin and cytochrome c 6 possess similar physicochemical and surface structural properties, which make them capable of replacing each other (Navarro, Hervas, & De la Rosa, 1997).
The monoclonal antibody Q425 is generally used to isolate CD4 receptors using immunoprecipitation techniques (Lynch, Dearden, Sloane, Humphery-Smith, & Cunningham, 1996). However, in a recent study, authors were surprised to see the effect of ethylenediaminetetraacetic acid (EDTA) on CD4-Q425 binding (Zhou, Hamer, Hendrickson, Sattentau, & Kwong, 2005). As the receptor CD4 does not utilize any metal (Healey et al., 1990), it was found that the antibody Q425 requires a divalent metal ion preferably Ca 2+ for its activity (Zhou et al., 2005).
In this study, we have performed a total of 2.45 μs of MD simulation on seven proteins replacing their cognate (bound in the structure) metal ions with another six similar ones. The results suggest that some of the proteins considered here can accommodate similar metal ions in their active site in the place of their cognate metal ions.

Materials and methods
2.1. Data collection and system setup In total, seven protein structures were used to perform the MD simulations. The details of the structures used in this study are provided in Table 1. The protein structures bound with a particular metal ion were extracted from the MIPS database (Hemavathi et al., 2010). The structures for MD simulation were chosen with the criteria that the structure (1) is highly resolved, (2) is of small size, (3) contains a single metal ion in the active site, (4) contains minimal set of other ligands, and (5) is published.
In total, seven proteins, as listed in Table 1, were chosen for MD simulations. In each case, the cognate metal (bound in the structure) was replaced with six other metal ions to study their binding to the protein. Thus, in total, 49 MD simulations were performed to study the binding of metals to different proteins. All crystallographic water molecules except those coordinating with active site metal ion and small molecules bound, if any, were removed from the protein models before MD simulations.

MD simulation
MD simulations were performed using the program GROMACS v.4.5.5 (Pronk et al., 2013). The optimized potentials for liquid simulations-all atom (OPLS-AA) force field was used for all of the simulations (Jorgensen & Tirado-Rives, 1988). A dodecahedron box was generated using the module editconf of GROMACS with the criterion that the minimum distance between the solute and the edge of the box was at least 1.0 nm. The protein models were solvated with the simple point charge (SPC) water model using the program genbox available in the GROMACS suite. Chloride and sodium ions, wherever required, were used to neutralize the overall charge of the system. Periodic boundary condition was applied along the three spatial directions to prevent edge effects. Energy minimization was performed using the steepest descent method with a maximum force cutoff of 1000 kJ/mol/nm. The solvent and ions around the protein were equilibrated in two phases. The first phase was performed under an NVT (canonical or isothermal-isochoric) ensemble for 100 ps with a reference temperature of 300 K. The second phase of equilibration was performed under an NPT (isothermal-isobaric) ensemble for 100 ps with a reference pressure of 1 bar. The reference temperature (300 K) was controlled by using a velocity rescaling thermostat (Bussi, Donadio, & Parrinello, 2007) with coupling constant .1 ps. The Parrinello-Rahman barostat (Parrinello & Rahman, 1981) was used to control the reference pressure 100 kPa with coupling constant 2 ps. The particle mesh Ewald method was used to compute long-range electrostatic interactions (Darden, York, & Pedersen, 1993;Essmann et al., 1995). A cutoff of 1.0 nm was used to compute the short-range van der Waals interactions. Bond lengths were constrained with the P-LINCS algorithm (Hess, 2008), and a time step of 2 fs was used to integrate the equations of motion. MD was performed for a time period of 50 ns for all 49 simulations discussed in this study. The protein-metal interaction energy was calculated using the equation where E protein-metal denotes the interaction energy between protein and metal. 'elec' and 'vdw' denote the electrostatic and van der Waals components of the energy, respectively.

Data analysis
Figures were generated using the program PyMOL (DeLano Scientific). Electrostatic potentials were calculated using the APBS (Baker, Sept, Joseph, Holst, & McCammon, 2001) module plugged into PyMOL. Most of the MD analyses were performed using the programs available in GROMACS and locally developed shell scripts. Graphs were prepared using the program Xmgrace (Turner, 2005).

Results and discussion
3.1. Bos taurus PLA 2 shows binding to most of the metals except Ba 2+ and Sr 2+ The cognate metal ion Ca 2+ of BtPLA 2 is heptacoordinated by carbonyl backbone oxygen atoms of Tyr28, Gly30, and Gly32, O δ1 and O δ2 of Asp49, and two water molecules ( Figure 1(A)). To study the replaceability of the cognate metal ion Ca 2+ of BtPLA 2 , we performed MD simulations of the protein by substituting its cognate metal with six other divalent metal ions such as Mg 2+ , Zn 2+ , Fe 2+ , Cu 2+ , Ba 2+ , and Sr 2+ . A root-meansquare deviation (RMSD) analysis shows that all the structures including mutants are well equilibrated during MD simulation ( Figure S1(A)). It has been shown that Ca 2+ is an essential cofactor for the functionality of the enzyme PLA 2 (Davidson & Dennis, 1990); in the absence of Ca 2+ , the enzyme loses its complete enzymatic activity (Dijkstra et al., 1981). Thus, it was tempting to understand whether the divalent Ca 2+ can be substituted by other similar divalent metal ions. An analysis of the interaction energy between the active site residue Asp49 and the bound metal ion suggest that Mg 2+ is the most preferred metal in the active site of the BtPLA 2 with an estimated binding energy of −625 kJ/mol (Figure 1(B)). The binding energy of the cognate Ca 2+ ion varies in the range from −380 kJ/mol to −250 kJ/mol. The metal ion Zn 2+ is stabilized in the active site with an energy of −360 kJ/mol, while the metal ions Fe 2+ and Cu 2+ remained bound to the active site of the enzyme throughout the simulation with energy of −350 and −450 kJ/mol, respectively. On the other hand, the analysis suggest that two metal ions Ba 2+ and Sr 2+ shift from the active site immediately within 2-3 ns of the MD simulation to the bulk solvent and thus show no binding to the BtPLA2 (Figure 1(B)). This might be due to the fact that Ba 2+ and Sr 2+ have larger metal ionic radius compared to other metal ions and thus cannot be accommodated in the active site. To further understand the dynamics of metal ions in the active site of the protein, we calculated the distance between the residue Asp49 and the metal ion during the MD simulation. The analysis shows that all the metals, except Ba 2+ and Sr 2+ , remain bound in the active site with an average distance of 2.5 Å. The metals Ba 2+ and Sr 2+ , however, move away from the active site and become part of the solvent phase ( Figure 1(C)).
To further understand the dynamics of active site geometry of the enzyme upon replacing the cognate metal ion with other divalent metal ions, we calculated the distance between Asp49 and all other metal-coordinating residues throughout the simulation. The active site residues show high fluctuation in the case of Ca 2+ -, Fe 2+ -, and Cu 2+ -bound protein, whereas, the protein bound to Mg 2+ , Zn 2+ , Ba 2+ , and Sr 2+ show low fluctuation (Figure 1(D)). Interestingly, in the case of Ba 2+ -and Sr 2+ -substituted structures, the active site residues remain rigid even in the absence of the metals. This suggests that the active site geometry of the enzyme BtPLA 2 is preformed to bind the divalent metal ions. To further investigate the residue-wise fluctuation of the protein molecule during the MD simulations, we computed the average root-mean-square fluctuation (RMSF) of all the seven structures of BtPLA 2 . In most of the cases, five set of residues (15-22, 30-35, 60-70, 77-87, and 107-123) exhibit high fluctuation with an RMSF of more than 1.0 Å ( Figure S1(B)). Interestingly, all these regions either belong to loops or highly solvent-accessible parts of the protein. Thus, the substitution of the cognate metal ion Ca 2+ by divalent ions shows negligible effect on the protein dynamics as a whole.

The enzyme SPP can accommodate most of the divalent metal ions
As reported in the literature, the enzyme SPP requires Mg 2+ ion for its activity. To further explore the possibility of the replacement of the cognate metal ion in the active site of the enzyme, seven MD simulations each of 50 ns were performed by substituting Mg 2+ with Ca 2+ , Zn 2+ , Fe 2+ , Cu 2+ , Ba 2+ , and Sr 2+ . The metal ion Mg 2+ is hexa-coordinated by three protein residues (two carboxyl atoms of Asp9 and Asp186 and one carbonyl oxygen atom of Asp11) and three water molecules (Figure 2(A)). An analysis of RMSD values of all the seven structures suggests that MD simulations were well equilibrated ( Figure S2(A)). Analysis of distance and interaction energy between the metal ion in the active site of the enzyme and its coordinating residue Asp9 reveals that all the substituted metal ions remain in the active site during the simulation (Figure 2(B) and (C)). Although, two metal ions (Ba 2+ and Sr 2+ ) show low binding energy, all other metal ions remain bound in the active site and exhibit high interaction energy (Figure 2(C)) suggesting that the enzyme can accommodate most of the divalent metal ions considered in this study.
It is observed that in the empty (without any other ligand e.g. sucrose or glucose) closed conformation of SPP, the cognate metal ion Mg 2+ is bound in the socalled non-catalytic site (Figure 2(A)). In this position, Mg 2+ is bound to Asp9, Asp186, Ser187, Asn189, and Asp190, but no longer with Asp11. To investigate whether the metal ions move toward the non-catalytic site during MD simulations, we calculated the distance between the carbonyl oxygen atom of the residue Asp11 and the metal ion. All the metal ions, except Zn 2+ , remain bound to the catalytic site of the protein molecule (Figure 2(D)). The metal ion Zn 2+ moves away from the catalytic site after~30 ns of MD simulation. When we further investigated whether Zn 2+ has moved to the noncatalytic site, we were not able to locate the metal ion near this position. It is to be noted that Mg 2+ binds to the non-catalytic site only in the closed conformation of the protein. Thus, the binding of the metal ions in the The active site geometry of BtPLA 2 . All the coordinating amino acid residues are shown as ball-and-stick model. The residue Asp49 essential for stabilizing the metal ion in the active site is colored in green. The cognate metal ion Ca 2+ and two coordinating water molecules are shown as sphere. In the mutant structures, Ca 2+ is replaced by Mg 2+ , Zn 2+ , Fe 2+ , Cu 2+ , Ba 2+ , and Sr 2+ in silico. (B) The interaction energy between the active site residue Asp49 and the metal ion bound. The interaction energy during the MD simulation for the cognate metal ion Ca 2+ is shown in black. The interaction energy for the substituted metal ions is shown in different colors as provided in the inset picture. The analysis reveals that Mg 2+ binds to the enzyme with high affinity than Ca 2+ itself. On the other hand, Ba 2+ and Sr 2+ show no binding toward the enzyme BtPLA 2 . (C) The distance between the active site residue Asp49 and the metal ion during the MD simulation. The two metals Ba 2+ and Sr 2+ which falls apart from the active site of the protein are shown in magenta and orange color, respectively. All other divalent metal ions remain bound in the active site with an average distance of 2.5 Å. (D) The dynamics of active site residues of the enzyme BtPLA 2 throughout the MD simulations. The graph suggests that even in the absence of metal ion, the active site geometry remains intact suggesting the pre-formation of the metal binding site. catalytic site is most likely due to the open conformation of the protein molecule.
An analysis of RMSF of the structures during MD simulations suggests that the residues involved in metal ion coordination in the active site of the protein molecule remain stable. However, regions (125-135, 150-153 and 215-222) of the protein molecule show high fluctuation during the simulation ( Figure S2(B)). These regions belong to the cap and the core domains, which are involved in the opening and closing process, of the protein molecule. It is shown that in the closed conformation, the cap domain moves toward the core domain by almost 25° (Fieulaine, Lunn, Borel, & Ferrer, 2005). It is to be noted that the protein transits from open to close conformation. However, it is not clear whether the process of opening and closing of the protein molecule is substrate dependent. Interestingly, an analysis of the dynamics of the active site of seven structures during MD simulations suggest that the substitution of the cognate metal ion (Mg 2+ ) with other similar divalent metal ions has no adverse effect on the geometry of the active site of the enzyme ( Figure S2(C)). The active site geometry of SPP. All the coordinate amino acids are shown as ball-and-stick model. The residue Asp9 essential for stabilizing the metal ion in the active site is colored in green. The cognate metal ion Mg 2+ and three coordinating water molecules are shown as sphere. The water molecule shown as sphere in magenta color corresponds to the position of Mg 2+ binding in the so-called non-catalytic site. The amino acids colored in yellow provide the coordination to the Mg 2+ bound in the non-catalytic site of the protein molecule. In the mutant structures, Mg 2+ is replaced by Ca 2+ , Zn 2+ , Fe 2+ , Cu 2+ , Ba 2+ , and Sr 2+ in silico. (B) The interaction energy between the active site residue Asp9 and the metal ion bound. The interaction energy during the MD simulation for the cognate metal ion Mg 2+ is shown in black. The interaction energy for the substituted metal ions is shown in different colors as provided in the inset picture. The analysis reveals that almost all the metal ions bind to the enzyme with high affinity. On the other hand, Ba 2+ and Sr 2+ show relatively lower binding affinity toward the enzyme SPP. (C) The distance between the active site residue Asp9 and the metal during the MD simulation. Almost all the metal ions remain bound in the active site of the enzyme during MD simulations. (D) The distance between carbonyl oxygen atom of Asp11 and metal bound in the catalytic site of the protein molecule. All metal ions, except Zn 2+ , are bound in the catalytic site of the enzyme. The analysis suggests that Zn 2+ moves away from catalytic site after 30 ns of MD simulation. However, the metal ion does not seem to occupy the non-catalytic site.

Pyrazinamidase binds almost all divalent metal ions
The Zn 2+ ion in the active site of PZase enzyme is hexacoordinated by O δ2 of Asp52, N ε2 of His54, and N ε2 of His71, and three water molecules in an octahedral geometry (Figure 3(A)). It is interesting to note that Zn 2+ ions mostly assume tetrahedral coordination liganded by sulfur atoms. However, octahedral coordination of Zn 2+ is also observed in the case of oxygen and nitrogen coordination (Bentley, Dodson, Dodson, Hodgkin, & Mercola, 1976). To investigate the binding of Zn 2+ ion in the active site of the PZase enzyme, we performed MD simulation for 50 ns of the structure bound with the cognate ion Zn 2+ . To further understand the substitutability of the cognate Zn 2+ by other divalent metal ions, we replaced Zn 2+ with Ca 2+ , Mg 2+ , Fe 2+ , Cu 2+ , Ba 2+ , and Sr 2+ and performed MD simulations for 50 ns each. The RMSD graph of all these seven simulations including the wild-type structure has been shown in Figure S3(A). The graph suggests that all the simulations were well equilibrated.
To see whether the cognate metal ion Zn 2+ and the substituted metals remain bound in the active site of the protein molecule during MD simulations, we calculated the interaction energy and the distance between the active site residue Asp52 and the metal ion. The analysis Figure 3. (A) The active site geometry of PhPZase. All the metal-coordinating amino acids are shown as ball-and-stick model. The residue Asp52 essential for stabilizing the metal ion in the active site is colored in green. The cognate metal ion Zn 2+ and three coordinating water molecules are shown as sphere. In the mutant structures, Zn 2+ is replaced by Ca 2+ , Mg 2+ , Fe 2+ , Cu 2+ , Ba 2+ , and Sr 2+ in silico. (B) The interaction energy between the active site residue Asp52 and the bound metal ion. The interaction energy during the MD simulation for the cognate metal ion Zn 2+ is shown in black. The interaction energy for the substituted metal ions is shown in different colors as provided in the inset picture. The analysis reveals that Mg 2+ binds to the enzyme with high affinity than Zn 2+ itself. On the other hand, Ba 2+ and Sr 2+ show lower binding affinity toward the enzyme PhPZase. (C) The distance between the active site residue Asp52 and the metal during the MD simulation. All the metal ions remain bound to the protein during MD simulations. All other divalent metal ions remain bound in the active site with an average distance of 2.5 Å. (D) The dynamics of active site residues of the enzyme PhPZase throughout the MD simulations. The graph suggests that the active site geometry remains intact during MD simulation. of interaction energy exhibits that most of the ions except Ba 2+ and Sr 2+ have high binding affinity. Ba 2+ and Sr 2+ bind the protein with lower affinity in the range from −300 to −200 kJ/mol. Interestingly, even though Ba 2+ and Sr 2+ have low binding affinity, they remain bound in the active site of the protein molecule (Figure 3(B)). Surprisingly, it is observed that Mg 2+ has better binding affinity than that of Zn 2+ in the active site of the PZase enzyme. This may be due to the fact that Mg 2+ is more stabilized by oxygen and nitrogen atoms. Analysis of distances between the active site residue Asp52 and the bound metal ion using the MD trajectories of all the wild-type and mutant structures suggests that all the divalent metal ions remain bound in the active site during MD simulations (Figure 3(C)). As it has been described that the binding of metal ion exhibits no perceptible conformational change in the protein molecule, the substitution of Zn 2+ with other divalent metal ion does not seem to affect the binding. Furthermore, as the metal ion is supposed to provide nucleophilic water for catalysis, any divalent metal ion may play the role of general base in this case. It is interesting to note that though the active site residues coordinating the metal ions belong to a loop of the protein molecule, the active site exhibits high rigidity during MD simulations (Figure 3(D)). It has been reported that the binding of Zn 2+ shows no detectable conformational change of the protein molecule, even in the metal-coordinating residues (Du et al., 2001). Although, all the loops (10-23, 52-72, 94-103 and 128-133) surrounding the active site of PZase show high flexibility during simulation ( Figure S3(B)), the metal-coordinating residues remain rigid. Surprisingly, in the cases of Fe 2+ -, Sr 2+ -, and Cu 2+ -bound structures, the regions 95-110, 105-115, and 155-165, respectively, show high flexibility; the reason, however, for this is not clear.

CDO does not bind to any divalent metal ions
To study the binding of the reported metal ion iron in its ferrous form (i.e., Fe 2+ ) and its substitution with other divalent metal ions to the enzyme CDO, we performed multiple MD simulations each of 50 ns using the program GROMACS. The proposed metal ion Fe 2+ is tetracoordinated by three N ε2 atoms of the residues His86, His88, and His140 along with a water molecule in an approximate tetragonal geometry (Figure 4(A)). The RMSD graph depicts that all the structures including mutants are well equilibrated during MD simulation ( Figure S4(A)). To see whether Fe 2+ binds to the enzyme, we calculated the interaction energy between the metal ion and the residue His140 during the MD simulation and found that the metal ion does not bind (Figure 4(B)). We also attempted to see if the metal ion remains bound to the protein molecule, we calculated the distance between Fe 2+ and active site residue His140 and observed that the metal ion shifts from the protein molecule to the bulk solvent with an average distance of 20 Å (Figure 4(C)). To further inquire whether any other divalent metal ion can be substituted in the place of Fe 2+ , it was replaced with Ca 2+ , Mg 2+ , Zn 2+ , Cu 2+ , Ba 2+ , and Sr 2+ and performed the simulation each of 50 ns again. The analysis of interaction energy and the distance between the metal and the residue His140 shows that none of these metal ions bind to the protein molecule (Figure 4(B) and (C)). This suggests that metal ion might not be required for the activity of CDO. However, the binding of the metal ion in the active site of the protein in the presence of the substrate or at a particular pH cannot be ignored either. It is to be noted that of three histidine residues, only His88 is protonated at N ε2 , whereas the residues His86 and His140 are deprotonated at N ε2 during the MD simulations.
In expectation that the dynamics of the active site geometry will fluctuate significantly during MD simulation, we calculated the distance between two groups containing His140 as one group and His86 and His88 as another group. We find that the fluctuation between these two groups varies between 2 and 3 Å, thus, suggesting a rigid active site even in the absence of the metal ion (Figure 4(D)). To further investigate whether the absence of metal ion in the active site has any adverse effect on the coordinating residues, we calculated the RMSF of each amino acid and plotted with respect to time of MD simulation. The regions 1-5, 20-25, 38-42, 64-66, 105-120, and 185-190 were found to show high fluctuation ( Figure S4(B)). All the three residues (His86, His88, and His140) supposed to coordinate with the metal ion do not exhibit much fluctuation. Thus, it can be suggested that CDO might not require a metal ion for its functionality.

Plastocyanin does not bind any divalent metal ion
Plastocyanin is a Cu 2+ -binding protein involved in electron-transfer shuttle from cytochrome b 6 f to P700 + of photosystem I in chloroplasts of algae, higher plants, and also in many cyanobacteria (Freeman & Guss, 2001). In the case of algae and cyanobacteria lacking plastocyanin proteins, cytochrome c 6 performs this function (Sandmann et al., 1983). The metal ion Cu 2+ is coordinated by four protein residues (His37, Cys84, His87, and Met92) in an irregular (or distorted) tetrahedron ( Figure 5(A)). To make sure that all the seven MD simulations were properly run, we calculated the RMSD values comparing the energyminimized crystal structure. The RMSD graph suggests that all the simulations were equilibrated ( Figure S5(A)). To calculate the protein-ligand interaction, the active site residue Cys84 was chosen and the interaction energy between the residue and the metal ion during the MD simulations was computed. The analysis exhibits that none of the metals considered in this study bind to the protein during MD simulations ( Figure 5(B)). A similar pattern can be found in the protein-ligand distance calculation ( Figure 5(C)). To further explore whether the active site of the protein changes due to the absence of metal, we computed the active site fluctuation during MD simulation. The analysis shows that the active site geometry of the protein remains rigid during simulations suggesting no structural perturbation due to the loss of the metal (Figure 5(D)). To further inquire whether the absence of metal ion in the active site has any adverse effect on the metal-coordinating residues, we calculated the RMSF of each amino acid and plotted with respect to time of MD simulation. The graph reveals that most of the residues show insignificant fluctuation with an average RMSF of 1.0 Å ( Figure S5(B)). It suggests that even in the absence of metal ion in the active site, the overall tertiary structure of the protein remains intact.

Anti-CD4 antibody Q425 binds most of the divalent metals except barium
Although one-third of all proteins are estimated to be metalloproteins, not many examples are known where antibodies interact with their cognate receptors with the help of metal ions. A recent study reported that the binding affinity of CD4 to its antibody Q425 increases in the presence of divalent metal ions at the interface The active site geometry of CDO. All the metal-coordinating amino acids are shown as ball-and-stick model. The residue His140 essential for stabilizing the metal ion in the active site is colored in green. The cognate metal ion Fe 2+ and the coordinating water molecules are shown as sphere. In the mutant structures, Fe 2+ is replaced by Ca 2+ , Mg 2+ , Zn 2+ , Cu 2+ , Ba 2+ , and Sr 2+ in silico. (B) The interaction energy between the active site residue His140 and the bound metal ion. The interaction energy during the MD simulation for the cognate metal ion Fe 2+ is shown in black. The interaction energy for the substituted metal ions is shown in different colors as provided in the inset picture. The analysis reveals that none of the metals bind to the enzyme. (C) The distance between the active site residue His140 and the metal during the MD simulation. None of the metal ion remains bound to the protein during MD simulations. (D) The dynamics of active site residues of the enzyme CDO throughout the MD simulations. The graph suggests that even in the absence of metal ion in the active site of the enzyme, the active site geometry remains intact suggesting the pre-formation of the metal binding site. (Zhou et al., 2005). The binding affinity between CD4 and Q425 increases by 55,000-fold higher in the presence of Ca 2+ compared with that of no metal ions. In this study, we performed the simulation of anti-CD4 antibody Q425 with bound Ba 2+ , and subsequently, Ba 2+ was replaced with Ca 2+ , Mg 2+ , Zn 2+ , Fe 2+ , Cu 2+ , and Sr 2+ in order to understand their relative binding to the protein. The cognate metal ion Ba 2+ is hexa-coordinated by three protein atoms (Asp32, Glu50, and Ser91) and three water molecules (Figure 6(A)). The overall quality of the MD simulation seems to be equilibrated ( Figure S6 (A)). The interaction energy between the bound metal and the active site residue Glu50 shows that Ba 2+ does not exhibit significant affinity toward the protein compared to other divalent metal ions (Figure 6(B)). As also observed experimentally, the metal ion Ca 2+ was found to bind with higher affinity (−370 kJ/mol) than Ba 2+ (−20 kJ/mol) (Zhou et al., 2005). Interestingly, the metal ion (Ba 2+ ) bound in the crystal structure has high temperature factor (63.53 Å 2 ) which is generally found to be in the range of 20-40 Å 2 . Although no biochemical study with Cu 2+ was performed, the metal ion shows the highest affinity (−430 kJ/mol) toward Q425 according to this study. In summary, Ba 2+ can be replaced by other divalent metal ions such Cu 2+ , Mg 2+ , Ca 2+ , Zn 2+ , Fe 2+ , and Sr 2+ in decreasing order of preference (Figure 6(B)).
To our expectation, the relative binding affinities of these metals match with the findings obtained in the Figure 5. (A) The active site geometry of plastocyanin. All the metal-coordinating amino acids are shown as ball-and-stick model. The residue Cys84 essential for stabilizing the metal ion in the active site is colored in green. The cognate metal ion Cu 2+ is shown as sphere. In the mutant structures, Cu 2+ is replaced by Ca 2+ , Mg 2+ , Zn 2+ , Fe 2+ , Ba 2+ , and Sr 2+ in silico. (B) The interaction energy between the active site residue Cys84 and the bound metal ion. The interaction energy during the MD simulation for the cognate metal ion Cu 2+ is shown in black. The interaction energy for the substituted metal ions is shown in different colors as provided in the inset picture. The analysis reveals that none of the metals binds to the enzyme. (C) The distance between the active site residue Cys84 and the metal during the MD simulation. None of the metal ion remains bound to the protein during MD simulations. (D) The dynamics of active site residues of the enzyme plastocyanin throughout the MD simulations. The graph suggests that even in the absence of metal ion in the active site of the enzyme, the active site geometry remains intact suggesting the pre-formation of the metal binding site.
original study (Zhou et al., 2005) which showed the highest affinity toward Ca 2+ (1.0), Sr 2+ (.58), and Cd 2+ (.24). All other divalent metal ions Mg 2+ , Ba 2+ , Mn 2+ , Co 2+ , Ni 2+ , Cu 2+ , and Zn 2+ showed lower or nondetectable affinities (<.005). However, it is to be noted that the binding affinity study was performed in the presence of 1 mM EDTA which is sufficient enough to chelate any weak-bound metal ions (Zhou et al., 2005). Thus, we propose that except Ba 2+ , all other divalent metal ions, which we considered in our study, can bind to the protein. A similar property of these metal ions affinity toward the protein is reflected in the calculation of protein-metal distance (Figure 6(C)); the figure shows that except Ba 2+ , all other divalent metal ions remain bound to the active site of the protein during MD simulation with an average distance of 2.3 Å from the active site residue Glu50 (Figure 6(C)). To further investigate the structural perturbation of the active site residues, we calculated the distance between two groups of the active site residues throughout MD simulation. The graph shows that structure bound with Ba 2+ exhibits high fluctuation (Figure 6(D)); this is most likely due to the detachment of the metal ion during MD simulation. Furthermore, the structures bound to Sr 2+ , Cu 2+ , and Mg 2+ show high variation toward the end of the simulation (Figure 6(D)). However, even with higher fluctuation in the active site geometry, metal ions remain bound to the protein. Thus, the flexibility of the active site residue Figure 6. (A) The active site geometry of anti-CD4 antibody Q425. All the metal-coordinating amino acids are shown as ball-andstick model. The residue Glu50 essential for stabilizing the metal ion in the active site is colored in green. The cognate metal ion Ba 2+ is shown as sphere. In the mutant structures, Ba 2+ is replaced by Ca 2+ , Mg 2+ , Zn 2+ , Fe 2+ , Cu 2+ , and Sr 2+ in silico. (B) The interaction energy between the active site residue Glu50 and the bound metal ion. The interaction energy during the MD simulation for the cognate metal ion Ba 2+ is shown in black. The interaction energy for the substituted metal ions is shown in different colors as provided in the inset picture. The analysis reveals that all divalent metal ions except Ba 2+ bind to the enzyme. (C) The distance between the active site residue Glu50 and the metal during the MD simulation. All the metal ions except Ba 2+ remain bound in the active site of the protein during MD simulations. (D) The dynamics of active site residues of the protein anti-CD4 antibody Q425 throughout the MD simulations. The graph suggests that the active site topology is highly flexible during MD simulation even in the presence of metal ions. It is most likely due to the absence of its binding partner, CD4 receptor. seems to be the intrinsic property of this protein. To further explore the flexibility of the residues, we averaged out the fluctuations of residues during the MD simulation and plotted them. The graph suggests that the regions 1-5, 40-60, 95-100, and 120-130 are highly flexible ( Figure S6(B)). Interestingly, the temperature factor for these regions in the crystal structure is higher than that of the average (48.9 Å 2 ).

Synaptotagmin 1 binds almost all the divalent metal ions
Synaptotagmin (Syt), a membrane-trafficking protein, contains an N-terminal transmembrane domain and two tandem cytoplasmic domains (C 2 A and C 2 B). Of 15 members of this family, eight  bind Ca 2+ with different affinities (Yoshihara et al., 2003) and are known to act as Ca 2+ sensors in synaptic vesicle docking and synaptic vesicle fusion (Fukuda et al., 2000;Pang et al., 2006). Sr 2+ , an analog of Ca 2+ , can also trigger neurotransmitter release, though less efficiently (Shin et al., 2003).
The Sr 2+ ion is hexa-coordinated by four protein residues (Asp303, Asp309, Asp363, and main chain oxygen of Tyr364), and one water molecule, thereby adopting an approximate tetragonal bipyramidal geometry (Figure 7(A)). The loop I which participates in coordinating metal ion is known to be destabilized in the presence of Mg 2+ ion whereas becomes stable in the presence of Ca 2+ or Sr 2+ Figure 7. (A) The active site geometry of synaptotagmin 1 C 2 B domain. All the metal-coordinating amino acids are shown as balland-stick model. The residue Asp363 essential for stabilizing the metal ion in the active site is colored in green. The cognate metal ion Sr 2+ is shown as sphere. In the mutant structures, Sr 2+ is replaced by Ca 2+ , Mg 2+ , Zn 2+ , Fe 2+ , Cu 2+ , and Ba 2+ in silico. (B) The interaction energy between the active site residue Asp363 and the bound metal ion. The interaction energy during the MD simulation for the cognate metal ion Sr 2+ is shown in black. The interaction energy for the substituted metal ions is shown in different colors as provided in the inset picture. The analysis reveals that all the metal ions bind to the enzyme. (C) The distance between the active site residue Asp363 and the metal during the MD simulation. All of the metal ions remain bound to the protein during MD simulations. (D) The dynamics of active site residues of the protein synaptotagmin 1 C 2 B domain throughout the MD simulations. The graph suggests that the active site topology is highly rigid during MD simulation. (Cheng et al., 2004). Interestingly, a similar observation was obtained in MD simulation as well. The loop I region shows high fluctuation in the presence of Mg 2+ compared to other divalent metal ions ( Figure S7(B)).
To find out whether the native divalent metal ion Sr 2+ can be replaced with other divalent metal ions, we performed MD simulations of Syt protein bound to Sr 2+ which was replaced by Ca 2+ , Mg 2+ , Zn 2+ , Fe 2+ , Cu 2+ , and Ba 2+ . The overall quality, as measured by the RMSD, of all the simulations is quite acceptable ( Figure S7(A)). The MD trajectories were analyzed and found that all the divalent metal ions considered in this study remain bound to the protein. The interaction energy reveals that the highest binding comes from Mg 2+ and the lowest with Ba 2+ (Figure 7(B)). Interestingly, it was observed in the earlier study that Mg 2+ can replace Ca 2+ , but with some structural changes at the active site (Cheng et al., 2004). Although no experimental details about Ba 2+ is available in the literature, our analysis suggests that most of these divalent metals can replace the cognate Sr 2+ ion in the active site of the synaptotagmin protein. Furthermore, we analyzed the movement of metal ions in the active site of the protein molecule during MD simulation. The metals stay intact in the active site with varying distance from the active site residue Asp363 (Figure 7(C)). For example, the lowest average distance between Asp363 and the bound metal is below 2 Å for Mg 2+ , whereas for Ba 2+ , it is approximately 2.7 Å (Figure 7(C)). It is to be noted that in the case of Mg 2+ -bound structure, the metal remains intact even when loop I shows high fluctuation during MD simulation ( Figure S7(B)); the loop I is involved in stabilizing the metal ions in the active site of the protein molecule. To further investigate the fluctuation of active site residues, we calculated the distance between two groups of the residues coordinating with the metal ion. We find, to our surprise, that the active site geometry of the protein molecule stands highly rigid during MD simulation (Figure 7(D)). The distance between two groups remains close to 1.33 Å throughout the MD simulation.

Conclusion
A total of 2.45 μs of MD simulation suggests that metal ions in the active site of certain proteins can be replaced with chemically similar ones. The enzyme PLA 2 which requires Ca 2+ for its catalytic activity shows better binding with Mg 2+ . In addition, metal ions Cu 2+ , Zn 2+ , and Fe 2+ also exhibit similar binding affinities to that of Ca 2+ , whereas, Ba 2+ and Sr 2+ do not display any binding to the enzyme. In the Ba 2+ and Sr 2+ substituents, interestingly, the active site geometry of enzyme remains rigid during MD simulation; most likely due to the absence of any metal in the active site. Although metal ions Cu 2+ and Fe 2+ show a similar affinity to PLA 2 to that of Ca 2+ , the active site geometry of these mutants show high fluctuation suggesting non-compatibility in the active site. In the case of the enzyme SPP, results suggest that none of the metal ions (Ca 2+ , Zn 2+ , Fe 2+ , Cu 2+ , Ba 2+ , and Sr 2+ ) considered in this study bind to the enzyme with better affinity than that of Mg 2+ . It is also reflected in the pattern of active site geometry movement during MD simulation; the mutant structures exhibit high fluctuation compared to the native (Mg 2+bound) protein. Furthermore, we also demonstrate that none of the metal ions fixed at the catalytic metal binding site move to the non-catalytic metal binding site during MD simulation. Although the metal Zn 2+ moves away from the catalytic site after~30 ns of MD simulation, it is not detected in the non-catalytic site. The enzyme pyrazinamidase depicts binding to Mg 2+ and Cu 2+ better than its cognate metal ion Zn 2+ . It is interesting to note that the active site topology of the enzyme shows no fluctuation during MD simulation only in the case of Cu 2+ -bound structure. Even in its cognate metal ion bound structure, fluctuation of metal-coordinating residues is high.
On the contrary, it is surprising to note that the enzymes CDO and plastocyanin do not even bind to its cognate metal ion Fe 2+ and Cu 2+ , respectively, during MD simulation. None of the metal ions substituted in the active site of these enzymes are retained during MD simulation. In the case of CDO, active site residues exhibit high fluctuation in the absence of metal ions. On the contrary, the active site topology of the protein plastocyanin remains intact even in the absence of metal ions. This suggests that the active sites of these two enzymes behave differently. Interestingly, in anti-CD4 antibody Q425, all the metal ions except its cognate metal Ba 2+ , substituted at the site of metal binding, remain intact during MD simulation. It is shown experimentally that Ca 2+ has more effect on its function than Ba 2+ . Although all the metal ions remain intact in the active site of the protein during MD simulation, the active site of the protein shows high fluctuation, most likely due to the absence of its binding partner CD4. For the protein synaptotagmin, almost all the metal ions exchanged in the place of its cognate metal ion show higher binding affinities than its cognate metal ion Sr 2+ . Moreover, the active site topology of the protein remains absolutely rigid throughout MD simulation.
Although, the results obtained in this study through MD simulations are promising, the artifacts known to be introduced in MD simulation algorithms cannot be ignored. The first error which could have been incorporated in this study is the fact that the exact force field parameters for several metals are not exactly derived. Secondly, it is a known fact that the binding of metal ions to proteins depends on its environmental structural pH; however, no measurement of these information was performed in this study. Thus, it is tempting to speculate that although some observations of this study might turn out to be correct, experimental evidences are required to support these results.