Noncovalent interactions between the second coordination sphere and the active site of [NiFeSe] hydrogenase

Abstract

QM/MM studies on seven truncated models of the oxidized as-isolated state of the [NiFeSe] Hases have been undertaken in order to find out the influence of the residues on the second coordination sphere on the active site. The major interactions of the second coordination sphere with that of the active site concerns hydrogen bonds between the CN ligand and Pro420 and Ala421 residues, weak S⋯N interactions between the Cys_(brdg) and the Sec with the His82 and Arg422 residues. Five types of weak noncovalent interactions between the ligands and the residues in the second coordination sphere are classified based on RGD and QTAIM methods. These residues have been found to affect the electronic structure of the active site as evidenced by the computed IR spectral features. The study reveals that the inclusion of all the residues in the second and the third coordination sphere of the enzyme may not be necessary for an accurate description of the active site features. The non-covalent interactions have been found to be stabilizing the Sec in the conformers.

---

Introduction

Photobiological hydrogen production with the aid of hydrogenases (Hases) has been a hot area of research. Hases are enzymes that catalyze the dihydrogen cleavage in sulphate reducing bacteria.¹ Studies on the electronic structure based on quantum mechanical (QM) methods have helped us to understand the pure organometallic nature of its active site.^2–6

Several studies on the cluster models of various states of the Hases are also available.^7–9 Amara et al. reported the study of [NiFe] Hases where the protein environment have been explicitly treated.¹⁰ The study reveals that for the Ni-A state of the [NiFe] Hase, the protein environment increased the Ni–Fe distance. The FeS clusters were also included and their role in the proton transport has been studied.¹⁰ Siegbahn and his group performed an exhaustive study of the mechanistic details of the formation of the states Ni-A and Ni-B of the [NiFe] Hase using a cluster model containing 120 atoms purely with the QM method.^11,12 Their results showed that the putative OOH⁻ species binds to Ni in an end-on fashion which was in fact not agreeing with the experimental findings of a bridged structure. They have also contributed to the understanding of the enzyme action through a QM/MM method and a combined QM/MM and cluster approach.^11,13,14 Jayapal et al. reported a comprehensive study on the Ni-A and Ni-B states of the [NiFe] Hases where the whole enzyme has been modeled through the QM/MM procedure and geometries close to the crystal structures were obtained. The study also found out that the QM/MM studies, rather than the isolated cluster calculations, correlate better with X-ray data and the structure of the unready Ni-A and Ni-B states were correctly predicted to have a bridging OOH⁻ and OH⁻ respectively.¹⁵ A similar study of the oxidized states have been reported by Söderhjelm and Ryde.¹⁶ Ryde has also brought out interesting guidelines on obtaining accurate energies and methods of truncation in proteins in general.¹⁷

Recently, Kampa et al. reported the study of the Ni-C state of the [NiFe] Hase which includes only the ligands in the first and the second coordination sphere upto a total of 165 atoms.¹⁸ The above work specifically addresses the impact of the His residue and this residue upon protonation influences the hydrogen bonding and the spin density at the bridging Cys (cysteines). There are many reports in support of the accuracy of the QM/MM models than the QM models.¹⁸ It has been reported that inclusion of the whole protein backbone does not yield very accurate results in a few cases and a full enzyme QM optimization may not yield accurate energies of the enzyme.^19–21 In the present work, QM/MM models of truncated structures of the Ni-SI state of the [NiFeSe] Hases from D. vulgaris Hildenborough (DvH) have been studied.

While there are several QM/MM studies on the [NiFe] Hases there are none on the [NiFeSe] Hases. Here the active site that performs the catalytic role in dihydrogen activation is treated with QM and the rest of the protein backbone that constitutes the second coordination sphere with MM. In any two layer calculations the two parts can be treated as two subsystems within the proteins.²²

The subgroup of Hase which are Sec (selenocysteine) containing, cytoplasmic and F420 nonreducing are called the Vhu Hases.^1,23 The residues Glu, Ala, Leu, Ser and Gly were found to be conserved in all the Vhu Hases.²⁴ But, aliphatic residues near the active site were found to be lesser in DvH (D. vulgaris Hildenborough) and DmB (Dm. baculatum) than in the other Hases of the Vhu family. Though the reason for this is evolutionary, the intriguing interactions that stabilize these residues in the second coordination sphere of the enzyme have to be understood. More importantly, the system addressed herein exhibits a peculiar Sec conformation, the reason for the existence of which may lie in the orientation of the protein skeleton at its proximity.^24,25 Therefore a thorough understanding of the nature and influence of the second coordination sphere on the active site is needed.

The pathway for the H₂ uptake in the DvH in the oxidized state has been already addressed through the Molecular Dynamics study (MD) and D₂/H⁺ exchange mechanism by De Lacey et al.²⁶ A cage effect was found to be modulating the enzyme catalysis, which is primarily due to the protein structure surrounding the active site. It was also proved that H₂ and D₂ access the active site through different pathways. MD simulations also show that the enzyme has the ability to catalyze double isotope exchange of D₂ with proton faster than the single exchange contrary to the [NiFe] Hases. These observations clearly reveals that the bimetallic core as well as the protein backbone are playing a vital role in the mechanism of action of this enzyme.²⁶

The present study aims (i) to understand the interaction between the protein backbone and the active site ligands in contact with it using a cost effective model for the Ni-SI state of the [NiFeSe] Hase of the DvH, through the RDG (Reduced Density Gradient) and QTAIM (quantum theory of atoms in molecules) methods, and (ii) to understand the interaction of the conformations of the selenocysteine ligands reported in the crystal structure.²⁵

Features of the protein backbone as observed in the crystal structure and the residues in the models reported here

The DvH is constituted by two polymer subunits A & B. Subunit A contains 317 amino acid residues and B the large subunit has 495 amino acids. The large subunit hosts the active site. Domain features of the enzyme were found to be typical of the uptake Hases. The active site contains Ni–Fe linked through two Cys residues (78B & 492B) which form the thiolate bridges (see Fig. 1). Ni is coordinated to a Csw75B which contains S in the dioxygenated form. It is also connected to the Ile74B & Gly76B. Such oxidized forms are observed in other Hases also.²⁴ For instance, both the terminal Cys ligands are oxidized in D. vulgaris Miyazaki. The Gly76B is in turn linked to the bridging S of Cys78. The N-terminal of the Cys78 is connected further to Pro79-Thr80-Ala81-His82. The His82 is proximal to the thiolate bridges at a distance of 4.08 Å from the S of Cys78. A Cys492 was connected to the additional S of the selenocysteines (Sec489) through a short loop containing two amino acids (Gly491 and Leu490). The amino acids Ala420-Pro421-Arg422-Lys423 were found to be closer to the CN groups in the active site. Hence, these four were also included in the MM optimization. The chain connecting the two terminal ligands (Csw & Sec) to the bridging cysteines are oriented in such a way that both the extrinsic S & Se of Sec489 are freely coordinated to the Ni. The electronic structure and the topological properties of the bidentate mode of Se–S coordination to Ni peculiar to the class of [NiFeSe] Hase has been investigated in our earlier work.²⁷ The above amino acids that constitute the second coordination sphere can be grouped into four categories based on their positions (i) those present above the Ni–Fe bridge (ii) the histidine group below the thiolate bridge, (iii) amino acids that flank the CO and CN ligands and the (iv) those that link the bridging thiolates with the terminal thiolate ligand of Ni. In order to understand their interaction with the active site, we have chosen the three models LARG, MDM and SML. In LARG all the groups from (i to iii) are present, in MDM groups listed in (iv) only are present and in the SML model only the aminoacids that noncovalent interaction with the ligands of the active site are present. HyL and HyP are models with only hydrophilic and hydrophobic residues in the second coordination sphere irrespective of their positions.

Fig. 1 QM/MM optimized geometry of the active site structure and the aminoacid groups in the second coordination sphere that engulf it in the truncated structure LARG obtained from Ni-OX state of the [NiFeSe] Hase as reported in crystal structure with PDB ID: 2WPN. Hydrogens are removed for clarity.

Computational details

Initial geometries were obtained from the X-ray structure of the enzyme DvH reported recently in the literature at a resolution of 2.04 Å (PDB ID: 2WPN).²⁵ Seven different models have been considered for the QM/MM studies (see ESI Fig. S1a–e†).

In all the models studied here the MM regions have been treated with universal force fields (UFF)^26,28 and QM optimization has been done at the BP86/TZVP level.^29,30 BP86 is used here since (i) earlier works in [NiFeSe] Hase³¹ and several other works in [NiFe] Hase^14,32 use BP86 due to the robustness of the function in handling transition metals and protein structures. (ii) For the sake of comparison and calibration with our earlier work on the QM studies of the active site of [NiFeSe] Hase.²⁷ As inclusion of long range electron correlation is important for the description of systems with noncovalent interactions^33,34 calculations were also carried out with the M06 functionals.³³

A few atoms have been restricted to positions as in X-ray data in accordance to the usual methods reported earlier.¹⁰ These include freely rotatable methyl groups in the amino acids namely C48 in Ala420, C26 & 27 in Ala81 and C7 & 8 in Leu49. The initial geometries for the MDM, SML, Hyl and HyP structures were obtained from the X-ray structure. QM/MM optimizations have been carried out after pre-treating the X-ray structure of the protein as follows. The whole enzyme as retrieved from PDB was first examined thoroughly and the surfactant molecules present along with the crystal structure as additives during crystallization were deleted. Then a centroid was fixed with respect to the Ni–Fe in the active site. All amino acid residues present beyond ∼7 Å from the centroid were carefully removed to arrive at the truncated structure. Then hydrogens were added to the structure. Firstly, only the added hydrogens were optimized at the QM level. This geometry was used as the initial structure for the ONIOM calculations. This pretreatment method was adopted based on earlier works on [NiFe] Hases.^13,35,36 All the short contacts between the active site and the second coordination sphere in the models studied here were analyzed using the package Mercury.³⁷ Cambridge Crystallographic Data Centre was used to identify the short contacts within the truncated structure.³⁷ For the study of hydrogen bonds the amino acid residues concerned along with the active site have been treated within the QM level. In each non bonded interaction artifacts due to truncation have been carefully identified and eliminated from the study. Short contacts observed in each of the models have been listed in Tables 2 and 3. A “Short” contact is one in which the internuclear distance is less than the sum of the van der Waals radii of the constituent atoms, and implies a noncovalent interaction.³⁸ The short contacts are further analyzed using the QTAIM and RDG methods at the BP86. For the QTAIM studies of the noncovalent interactions, calculations were carried out at the M06 level also.

The large model (LARG) consists of 200 atoms of which 47 atoms in the active site have been treated as part of the QM system and the remaining residues in the MM layer as implemented in ONIOM in G09W.³⁹ LARG contains nine amino acid residues that constitute the second coordination sphere, including those encapsulating the ligands to the Fe, and those that connect the bridging cysteines (Cys_(brdg)) and the terminal residues in the Ni in the MM layer. For the LARG model, calculations were also performed at the M06 level where C, H, O & N were treated with the 6-31g(d), Ni & Fe atoms with the 6-311G(3df) level and S and Se atoms with the 6-311g(dp) basis sets. In the MDM structures 137 atoms were present of which 42 atoms were treated at the QM level (see Table 1). MDM-Conf-1, MDM-Conf-2 & MDM-Conf-3 represent the three conformers of the Sec were modeled. The model SML contains the residues near the CN groups and the His ligand below the thiolate bridges. The hydrophilic and hydrophobic residues forming the second coordination sphere of the active site have been considered as two different models namely HyL and HyP (see Fig. 8). HyL contains completely the hydrophilic residues in the second coordination sphere irrespective of the positions and HyP contains residues that are hydrophobic in nature.

Table 1 Models studied and the residues in the proximity of Ni, Fe or Sec included in the model as part of the MM layer. Residue included in the calculation are indicated as (↙) excluded are indicated as (×)

	Group I near Ni	Group II His	Group III near Fe	Group IV near Sec	Other residues present
LARG	↙	↙	↙	↙	—
MDM	↙	↙	×	×	—
SML	×	↙	↙	↙	—
HyL	↙	↙	↙	×	Glu
HyP	×	×	×	×	Ala444, 493, 420 & 81, Leu425 & 490, Val77 & Ile74

---

Results and discussion

Geometric parameters of the LARG, MDM-Conf-1, SML, HyL and HyP have been compared with the crystal structure in Table 2. Results show that the Root Mean Squared Deviation (RMSD) in the bond length parameters between the five models lay within 0.043–0.08. Upon optimization the Ni–Fe bond lengths generally increase in almost all cases by ±0.2 Å. This increase is also implied in the bond angles and dihedral. Metal–sulphur bridges are described well in all the models. Ni–Se, Ni–S and S–Se bond distances are constant in all the models (2.4, 2.2 and 2.2 Å) respectively. The Ni–Se, S–Se and Ni–S_(brdg) bonds match well with those in the crystal structure.

Table 2 Selected bond lengths (Å) and the root mean squared deviations of the geometric parameters of the models studied from that of the crystal structure²⁴

	X-ray	LARG BP86	LARG M06	MDM-Conf-1	SML	HyP	HyL
NiFe	2.481	2.680	2.677	2.286	2.702	2.716	2.481
Ni–S(O2)	2.275	2.295	2.293	2.338	2.295	2.322	2.160
Ni–Se	2.338	2.439	2.429	2.408	2.435	2.464	2.445
Ni–S(Se)	2.128	2.193	2.195	2.230	2.211	2.216	2.185
S–Se	2.127	2.194	2.194	2.199	2.199	2.202	2.127
Se–C49	1.959	2.014	2.014	2.023	2.026	1.988	1.955
C–N d	1.184	1.179	1.178	1.178	1.188	1.179	1.184
C–N p	1.150	1.174	1.174	1.177	1.188	1.178	1.151
C–O	1.124	1.169	1.168	1.166	1.178	1.168	1.125
Fe–C d	1.847	1.884	1.887	1.887	1.874	1.872	1.848
Fe–C p	1.841	1.858	1.857	1.888	1.877	1.887	1.842
Fe–C(O)	1.804	1.715	1.715	1.710	1.712	1.713	1.803
Ni–S(brdg)d	2.318	2.284	2.283	2.275	2.249	2.242	2.287
Ni–S(brdg)p	2.266	2.286	2.284	2.233	2.239	2.268	2.160
Fe–S(brdg)d	2.271	2.239	2.236	2.289	2.284	2.290	2.287
Fe–S(brdg)p	2.288	2.265	2.265	2.274	2.274	2.304	2.236
RMSD		0.049	0.043	0.073	0.078	0.823	0.052

---

The proximal and distal C–N bond lengths are unequal in crystal structure, whereas they are equal in all the models except the HyL and LARG. This may be due to the presence of residues closer to the proximal CNs only in these two models. These residues seem to increase the electron density localized near the CN ligands as they behave as donors in hydrogen bonds and so through the noncovalent interactions they increase the Fe–C bond lengths slightly.

The effect was not brought out in the MDM where these residues are absent. In HyP residues present to the CN are purely aliphatic in nature and do not have any noncovalent interaction with the CN group. Therefore this effect was not brought out in HyP. Similarly the Ni–S_brdg (distal and proximal) distances are not equal in the crystal structure.

The mean squared deviations in the computed CO bond lengths are also closer to that of the X-ray data only in HyL. Fe–S (proximal and distal) bond lengths are nearly equal in the crystal structure but are described less accurately in the HyP & HyL. In the LARG model this is overcome. The RMSDs are lower in LARG and HyL. Therefore the presence of secondary ligands connecting the bridging and the terminal Cys as in MDM-Conf-1 and the inclusion of all the hydrophobic residues as in HyP may not be mandatory in modeling the active site features. In the MDM-Conf-1, as the Fe–C proximal and distal bond distances are similar unlike the LARG and SML, the presence of the residues near the CN ligands certainly influences the nature of the Fe centre. The dihedral plane containing the S–Se–C–C of the Sec also differs by 10° from the X-ray data in MDM-Conf-1. This is due to the absence of the Arg422 residue closer to the Sec in MDM-Conf-1. Therefore the Arg422 residue influences the position and orientation of the Sec and hence influences the conformations in the Ni-SI state of the [NiFeSe] Hase.

Hydrophilic and hydrophobic residues

Residues in the second coordination sphere both hydrophobic and hydrophilic were identified and their role in dictating the spatial arrangement of the bimetallic core has been analyzed in this section. The residues His82, Arg422 and Glu28 were the only hydrophilic ones within the region of study. Residues were identified to be hydrophobic or hydrophilic based on the scale of hydrophobicity calculated using the code WebLab Viewerlite 3.2 from Molecular Simulations, which assigns hydrophobicity based on the residues surrounding a particular amino acid in a protein. Gly76 which is the fourth amino acid in the hydrophilic residues is present at a distance of 3.78 Å from the Csw75. As this is present in the third coordination sphere it has not been included in the HyL model.

Glu28 in HyL is a conserved residue for all Vhu and Fru Hases. Hence their presence must be of significance to the enzyme action. In DvH the Glu28 is placed between the proximal bridging Cys and the Csw at a distance of 3.65 Å from Csw suitable for the proton transport to Csw. Among the hydrophobic residues, (Ala81, 420, 444 & 493, Leu425 & 490, and Ile74) Ala493 & Ala444 are positioned at least 4 Å away from the active site and posses no direct noncovalent interactions with those in the active site (see ESI Table ST1†). Topological properties reveal the absence of hydrogen bonds. Ala420 & 81 are closer to the CO ligand in the active site where, the free β methyl groups of Ala420 & 81 flank the CO. While the Ala420 shares a hydrogen bond with the free C-terminal of the Arg442, Ala81 and has a clear short contact with the Cys78. The hydrophilic residues present within the region of study namely Arg442, His82, Gly76 and Glu28 were found to influence the orientation and position of the side chains of Cys, Sec & Csw. The α-C of the Sec falls within the van der Waals radii of the O of the Glu28. Glu28 is placed in such a way that one of the C-terminal groups faces the Arg442 and Sec. Three interactions shorter than the van der Waals distances were observed with the distal and central N of the Arg442 and the O of the carbonyl ligand.

Conformations of the Sec

In DvH the Sec is known to be present in three conformations.²⁵ Reports show that in DvH where Sec replaces a Cys ligand of DmB, has increased pH activity due to the presence of selenium atom. The oxygen tolerance is also suspected to be due to the conformations (1 & 2) where the approach of the O₂ and CO₂ to the Ni is blocked by the Sec orientation.²⁶ Therefore to understand how these conformations could be realized in the active site in the presence of residues in the second coordination sphere, the conformations 2 & 3 (MDM-Conf-2 & 3 respectively) have been modeled from conformer 1 (MDM-Conf-1). The MDM-Conf-3 has an extrinsic sulphur (S²⁻) lesser than the other two conformers in the Sec ligand.

CO and CN⁻ stretching frequencies

The presence of the two non-proteinaceous ligands CO & CN in the active site has been of great advantage in understanding the details of the oxidation state, ligand positions and the overall chemistry around the metal centers.⁴⁰ The computed CN stretching modes were coupled with each other and hence one higher and one lower frequency has been observed for the two CN ligands corresponding to their symmetric and anti symmetric modes (Table 3). Comparison of the CN stretching frequencies shows that the anti symmetric stretching of LARG excellently agrees with the experiment. While the other models MDM and HyP in particular show lower absorption frequencies than expected.

Table 3 Characteristic IR stretching frequencies of the CO & CN ligands in the models studied against the experimental frequencies (cm⁻¹)

	QM/MM extrapolated energy (a.u)	νCO in cm^−1	νCN in cm^−1
LARG	−7908.3218	1946	2084, 2124
SML	−7789.0105	1945	2066, 2086
MDM-Conf-1	−7585.0254	1883	2008, 2018
HyP	−7376.9527	1879	1999, 2012
From Exp.^24	—	1939	2079, 2085, 2094

---

CO stretching modes (Table 3) were scaled by an additive correction of 28 cm⁻¹ as suggested for metal carbonyl complexes^41,42 and also used widely in literature.¹⁸ The carbonyl stretching frequency (ν_CO) value in HyP and MDM match well with the experimental value. In LARG the ν_CO values are in good agreement with the expected value of 1939 cm⁻¹. This is due to the presence of ligands near the CN group which increases the electron localization around CN and reduces the electron density at Fe and obviously a reduced Fe to CO back bonding. Therefore a higher ν_CO value has been observed in this model than the models without these residues. In the case of MDM where the side chains responsible for CN electron localization are absent, IR frequencies lesser than those observed for the X-ray structure. In summary, analysis of the computed IR data shows that the LARG model portrays the Fe coordination environment well as found in the crystal structure.

Noncovalent interactions between the first and second coordination sphere residues

Steiner has proposed a definition for hydrogen bonds (H-bonds) which is widely applicable to a variety of chemical entities and encompasses a gamut of intermolecular interactions ranging from solids to protein chains with strength varying from 0.4 to 40 kcal mol⁻¹. Accordingly, “AD–H⋯A interaction is called a hydrogen bond, if the following conditions are true. (1) It constitutes a local bond, and (2) D–H acts as a proton donor to A”.⁴³ Based on the range of observed hydrogen bond strengths of proton donors and proton acceptors, Steiner has categorized H-bonds and shown that H-bonds exist in a broad range with a continuum of strengths.⁴⁴ For practical purposes he classified them into weak, strong and intermediate type. It has been proposed that the physical basis for these three types of H-bonds are different and therefore their energies. Maximum H-bond strength is achieved when the D–H⋯A angle is 180° and is relaxed considerably in the case of moderate and weak bonds.

Bond length and bond angle parameters that satisfy the standard criteria are searched for at these donor–acceptor sites and the results are analyzed in detail. Tables 4 and 5 lists out the significant short contacts in LARG model. A search through the various short contacts showed that the CN ligands of the Fe establishes four H-bonds with the two Hs of Arg422 one from the Pro421(C–H65⋯N) and one from the H of Ala(C–H68⋯N). The two H-bonds with Arg are due to the N donors (C=N⋯H62 and C=N⋯H81). They are labeled as interactions 1–4 in Fig. 2. C–H68⋯N is from Pro421. These interactions have A⋯H distances shorter than the sum of their van der Waals radii and correspondingly lengthened D–H bonds (Table 5).

Table 4 Classification of H-bonds within active site and the residues in the second coordination sphere in 2WPN based on Steiner's analysis²⁴

Residue⋯ligand		van der Waals radii (Å)	Bond length (Å)	Bond angle CX⋯H	Dihedral angle CXHY	Type of H2 bond
Arg422⋯CN	N43–H62	2.75	2.6227	132.52	−144.8	Dispersive
Arg422⋯CN	N43–H65	2.75	2.4266	84.51	66.9	Dispersive
Arg422⋯CN	N43–H68	2.75	2.7088	116.31	−149.8	Dispersive
Arg422⋯CN	N43–H81	2.75	2.3901	136.17	−53.2	Strong electrostatic
Sec	S36⋯N28	3.35	3.2241	84.06	−72.4	—
Sec	Se35⋯H57	3.10	3.0218	119.12	102.8	—
Csw	H50⋯O4	2.72	2.2733	152.18	174	—

---

Table 5 Details of hydrogen bonds between the active site and the residues in the second coordination sphere and the corresponding A (acceptor)/D (donor)⋯H bonds

Residue⋯ligand	A⋯H–D	Shortening of A⋯H (Å)		Lengthening of H–D bond (Å)
		R(A–H) van der Waals distances	R(A⋯H) observed	R(H–D) expected	R(H–D) observed
Arg422⋯CN	N43⋯H62–N19	2.75	2.623	1.000	1.044
Arg422⋯CN	N43⋯H65–C11	2.75	2.427	1.070	1.070
Arg422⋯CN	N43⋯H68–C18	2.75	2.709	1.070	1.110
Arg422⋯CN	N43⋯H81–N28	2.75	2.390	1.000	1.001
Sec	S36⋯N28	3.35	3.374	—	—
Sec	Se35⋯H57–C25	3.10	3.022	1.070	1.109
Csw	H50⋯O4–S2	2.72	2.273	1.070	1.529

---

Fig. 2 (a) Hydrogen bonds between CN ligand and second coordination sphere residues. (b) Molecular graph showing the five hydrogen bond interactions between CN ligand and the residues in the second coordination sphere, where bond paths with (3, −1) critical points are represented in red color.

All these A⋯H interactions are shorter than 3 Å and have a θ(DHA) up to 136°. The N–H81⋯N, bonds has the maximum directionality among the four and therefore must be stronger than the other types of H-bonds. In these H-bonds 1 & 3, the D–H bond lengthens up to 0.003 and 0.004 Å. Therefore they may be classified as dispersive H-bonds as prescribed by Steiner.⁴³ It is well known that the influence of H-bonds can be realized on the acceptor (A) side, weakening the acceptor bond (A–C) and lowering the IR stretching frequency. This is clearly revealed in the increase of the CN (proximal) distances and a consequent shortening of the Fe–C bonds in LARG has been observed. Hence the presence of the H-bond has weakened the CN stretching frequencies in this model (see Table 4).

Arg442 also exhibits an interesting Se⋯H interaction with the Sec. Se⋯H(D) interactions are also known in literature to be equal in strength to that of the O⋯H and S⋯H interactions. The H(57)⋯Se interaction between the Sec and the Arg were also found to be shorter than their van der Waals distances. The stabilization gained due to this can be attributed to the stability of the conformation of Sec. The Sec orientation in the Conf-1 in the active site is also stabilized by another strong H-bond with the γ-H to the Se of Sec and the O (4) of Csw. This O⋯H bond falls well within the strong H-bond distances predicted by Steiner. A six membered ring motif containing the donor (C) acceptor (Se) atom could be spotted due to this H-bond (Fig. 3). The proton donor being a C atom, in the Se⋯H H-bond, A⋯H distances shorten twice as that of the X–H lengthening (0.08 & 0.04 Å respectively). The H-bond between the Csw and Sec shows higher directionality than the other H-bonds (θ(DHA) ≈ 152°). The His82 is well known in Hases to play a major role in proton transfer in the active site.²⁶ This residue has been reported to be a conserved feature in [NiFeSe] Hases also. The N43 of His is in short contacts with the CO ligated to Fe (3.195 Å) and the H153 of the His was found at a distance of 3.36 Å with the S6 of the bridging thiolate (distal).

Fig. 3 Molecular graphs revealing the six membered motif formed due to hydrogen bond between H of Sec and O of the sulfinate from Csw75 in Conf-I.

Two S–N short contacts were also observed in LARG. These short contacts have been further quantified using the reduced gradient electron density and electron density topology. A CCDC investigation of the crystal structure also shows an S⋯N short contact existing between the Arg in the second coordination sphere and the Sec in the first coordination sphere. S⋯N contacts are a recent entry to the long list of Noncovalent interactions known in literature. First report of S⋯N contacts was made by Rosenfield and coworkers.⁴⁵ Rosenfield' classification of the S⋯N short contacts shows that the mean value for the ZSN angle in type-II contacts (were Z = nucleophile) would be 30° (see Fig. 4). In the present case this angle is 21.5°. J. S. Murray established that such contacts do occur in crystal structures to stabilize them and are weaker than the halogen bonds.^46–48 The S⋯N interaction, Se⋯H bond and the S–O⋯H H-bond between Csw75 and Sec stabilize the Conf-1 in the bidentate mode of S–Se bonding to Ni. This may be the reason why the Conf-1 is the most stable conformation. In the MDM-Conf-3 modeled in the present work, the H⋯O–S H-bond shortens further to 2.12 Å but the θ(H⋯OS) remains to be 149° similar to the MDM-Conf-1. While in MDM-Conf-2 this contact is lost due to the Sec orientation.

Fig. 4 Schematic representation of the S⋯N interaction. (a) Rosenfields type-II S⋯X close contacts show that the mean value of the ZSN angle in contacts (with Z = nucleophile) would be 30° (b) SeSN angle in Conf-1 of Ni-SI state of [NiFeSe] Hases.

The tendency of S to bond with the nucleophiles like N and O has been attributed to the presence of strongly directional σ-holes.⁴⁹ Scheiner has investigated such interactions in isolated molecules and reported that S⋯N bonds like the P⋯N, As⋯N and Cl⋯N can be made stable and comparable in strength to that of H-bonds if electronegative atoms like halogens are substituted at the sulphur.^49–51 One of the three minima structures obtained by Scheiner had only the S⋯N bonds and neither H-bonds or halogen bonds. In these cases X–S⋯N angle was found to be around 170° and S⋯N length as short as 2.4–2.7 Å. In the case of the LARG model studied here, the N65 of the Arg lies at distance of 3.22 Å which is only 0.13 Å shorter than the van der Waals radii. Though the formation of such contacts and their stabilization energies cannot be computed in truncated models of large enzyme systems, RGD analysis and the topological features in the next section help identify the nature of these interactions. In the HyP model the hydrophobic residues Ala81, 420, 444 and Ile425 flank the carbonyl group and display weak noncovalent interaction through oxygen. The residues Ala81, Ala420 & Ala444 are positioned one below the carbonyl, one at the distal and the other at the proximal faces to the carbonyl ligand. These amino acids posses free methyl groups that orient towards the oxygen and hence create a hydrophobic environment around the CO.

Reduced density gradient (RDG) analysis

Reduced density gradient methods have been in use for at least four decades. RDG analysis proposed by Yang et al.⁵² is an index for identifying and analyzing the noncovalent interactions (NCI) based on the electron density and its reduced gradients. In NCI analysis, all the interactions with at least a specified fraction of the density from within a molecule are turned off, screening out only the intermolecular interactions. The resultant intermolecular interactions are depicted as surfaces of various color codes-red signifies strong repulsion, green stands for weak interactions and blue signifies strong attractions.

NCI analysis reveals a set of complex interactions between the ligand and the protein, which arise from a combination of specific atom–atom interactions (i.e., H-bonds) as well as broad surfaces indicative of stabilizing van der Waals interactions. NCI analysis clearly highlights how a ligand supports the geometry of the active site, and the many weak interactions that contribute to the interaction energy between the ligand and protein. In MDM-Conf-1, the region in the centre of the three membered ring containing S–Se–Ni shows a strong attractive blue isosurface indicating bonding interaction and a ring structure. The H attached to the γ-N of the Sec shows a weak attractive isosurface with the α-C of the Sec (see Fig. 6a). Sec also has a significant (S)O⋯H(C) interaction which is classified as a weak donor acceptor type in literature. In Conf-2 such interactions are absent both with γ-N and with α-C (see Fig. 5). While in Conf-2, the CN group (distal) involves in a weak noncovalent interaction with the α-C of Sec. The Se–S bonds due to their orientation away from the thiolate bridges, are parallel to the SO₂ group and hence possess a very weak attractive interaction with the electronegative oxygen and the Se–S. In the large active site where the amino acids which make the second coordination sphere like the His placed below the thiolate bridges is present and they have a significant interaction with the bridging S and the hydrogens of the α-C. The Arg group present above the S–Se plane also shares attractive interactions with the N–H atoms of the protonated amino acids.

Fig. 5 (a) & (b) NCI plot showing attractive surface between Sec and Csw, where isosurfaces correspond to an S value of 0.3 a.u and RGB scale ranges from −0.1 < ρ < 0.1 a.u.

Fig. 6 (a) 3D-NCI plot showing attractive surfaces indicating a hydrogen bond between Sec and Csw in MDM-Conf-1 (b) absence of attractive isosurface in Conf-II indicates lack of hydrogen bond between Sec and Csw. Isosurface correspond to an S value of 0.3 a.u and RGB scale ranges from −0.1 < ρ < 0.1 a.u.

Quantum theory of atoms in molecules (QTAIM) analysis

Further to understand the nature of these interactions, a QTAIM analysis was performed on the LARG model and for the sake of clarity and considering the computational cost, the LARG model and the residues in it have been discussed as three fragments on which topological properties have been analyzed. Molecular graph features four key amino acid residues in the proximity of the distal CN and they share five different interactions which are present as bond paths originating from the N43 of the CN (distal) to H (57, 62, 65, 68 & 81) of the amino acid side chains. These interactions are labeled as 1–5 respectively (Fig. 2). Topological properties at the BCPs corresponding to these interactions are presented in Tables 6 and 7. A negative Laplacian value, a positive H(r) and a V/G ratio of less than 1 shows that these interactions are of transient type i.e. they can also be called-both weakly covalent and weakly ionic. That is, from the magnitude and sign of the Laplacian values these bonds should be called ionic but can also be called covalent, based on the positive sign and absolute values of H(r). While interactions 1, 2 & 3 are very weak, interactions 4 is strong. Their ρ(r) values at the BCPs in interaction 4 is indicative of hydrogen bonding interaction.⁵³

Table 6 Topological properties at the BCPs of H-bonds at BP86/TZVP (in a.u)

		ρ(r)	L(r)	G(r)	H(r)	V(r)	G(r)/ρ	V(r)/G(r)
Ala⋯CN	N43–H62	0.0073	−0.0059	0.0050	0.0009	−0.0041	0.6833	−0.82
	N43–H65	0.0125	−0.0105	0.0090	0.0015	−0.0074	0.7160	−0.83
	N43–H68	0.0069	−0.0054	0.0046	0.0008	−0.0038	0.6617	−0.83
	N43–H81	0.0118	−0.0099	0.0088	0.0011	−0.0077	0.7456	−0.87
Arg⋯Sec	S36–N28	0.0119	−0.0083	0.0076	0.0007	−0.0069	0.6377	−0.91
	Se35–H57	0.0070	−0.0054	0.0042	0.0012	−0.0030	0.5968	−0.72
Sec⋯Csw	H50–O4	0.0153	−0.0117	0.0117	−0.0001	−0.0118	0.7648	−1.01
Sbrdg⋯His	S21–N16	0.0120	−0.0084	0.0076	0.0008	−0.0068	0.6294	−0.90
	S6–H38	0.0035	−0.0027	0.0020	0.0007	−0.0013	0.5632	−0.66
	N16–CO	0.0073	−0.0062	0.0052	0.0011	−0.0041	0.7094	−0.80

---

Table 7 Topological properties at the BCPs of H-bonds at M06/TZVP (in a.u)

		ρ(r)	L(r)	G(r)	H(r)	V(r)	G(r)/ρ	V(r)/G(r)
Ala⋯CN	N43–H62	0.0071	−0.0049	0.0042	0.0008	−0.0034	0.5915	−0.81
	N43–H65	0.0121	−0.0093	0.0080	0.0019	−0.0061	0.6612	−0.76
	N43–H68	0.0069	−0.0050	0.0042	0.0006	−0.0036	0.6101	−0.86
	N43–H81	0.0114	−0.0093	0.0076	0.0017	−0.0056	0.6661	−0.74
Arg⋯Sec	S36–N28	0.0112	−0.0083	0.0073	0.0010	−0.0062	−0.6524	−0.85
	Se35–H57	0.0071	−0.0052	0.0043	0.0009	−0.0033	−0.6056	−0.77
Sec⋯Csw	H50–O4	0.0153	−0.0125	0.0109	0.0016	−0.0093	0.7124	−0.85
Sbrdg⋯His	S21–N16	0.0112	−0.0083	0.0073	0.0010	−0.0063	0.6518	−0.86
	S6–H38	0.0033	−0.0027	0.0021	0.0006	−0.0015	0.6287	−0.71
	N16–CO	0.0071	−0.0063	0.0052	0.0011	0.0041	0.7324	−0.79

---

H57 of Arg present at a distance of 3.9 Å from the Sec shows two weak BPs with the Se and S. H57 forms a weak hydrogen bond with the Se atom. His interactions with the S_brdg are also found to be weak. Fig. 7 shows two bond paths that originate from N22 represented as 1 & 3. A bent path 3 indicates a weak interaction with the CO ligand to Fe. Bond path 2 is an S⋯N interaction.

Fig. 7 Molecular graph showing interactions between His and the S_brdg, where bond paths with (3, −1) critical points are present.

Fig. 8 (a) QM/MM optimized geometries of the HyP model and (b) HyL models at BP86/TZVP.

[NiFeSe] Hases of DvH vs. DmB

Comparison of the residues in the second coordination sphere of the Ni–SI and Ni–R of the [NiFeSe] Hases of DvH and DmB respectively shows that the following features in the second coordination sphere are similar to both. [NiFeSe] Hases have a basic residue above the active site closer to the acidic Sec and an acidic His below the bridging Cys and a hydrophilic Glu residue responsible for proton transfer near the Sec and Csw, aliphatic residues closer to the CN & CO ligands. These features are also common to DmB and provide the clue that the Ni may be the probable site for proton activation suitably surrounded by hydrophilic and acidic residues.

Conclusions

The recently reported oxidized as-isolated state of the [NiFeSe] Hases of the enzyme DvH have been investigated through QM/MM method. Seven different truncated models have been studied to understand the features of the active site of the enzyme. While the model which includes residues closer to the CN & CO ligands have been found to be important in describing the IR spectral features of the enzyme, the side chains that link the terminal and bridged thiolate but are not directly linked to the metal sites were not found to alter the geometric features of the enzyme. Hence a truncation after the primary coordination sphere is much valid in [NiFeSe] Hases as in other Hases. This may well be an opportunity to investigate the active site from a QM perspective. Various non covalent interactions that stabilize the three conformations of Sec have been addressed through the NCI and QTAIM analysis. Formation of a six membered H-bonded motif of the γ-H of the Sec with the O4 of Csw and a S⋯N short contact with the Arg has been identified to be the major stabilizing factor for the stability of Conf-1. This may also hold the key for the oxygen tolerance of the enzyme. The MDM-Conf-2 was found to be stabilized by 77.73 kcal mol⁻¹. Additionally, the Arg C-terminal was found to be in short contacts with the CN (distal) resulting in a reduced stretching frequency. From the study, it is evident that the acidic, basic and proton accepting residues are present at the distances favorable for interaction within the second coordination sphere. The noncovalent interactions within these residues surely are exploited by the active site to the full benefit of the enzyme action.

A comprehensive outlook of the active site has been developed through this work. In the light of its importance in the oxygen tolerance, the ambiguity in the crystal refinement studies and the long unanswered question of the role of Se, this work has modeled the structure of the active site. While many sophisticated methods are available, the choice of the QM/MM method on small truncated models does offer opportunity for further understanding of the enzyme at an affordable cost and at the same time without compromising accuracy.

Acknowledgements

We thank the reviewers for their critical comments and suggestions. SAV thanks the University Grants Commission (UGC), India and Bishop Heber College for FDP program (Ref. No. F.ETFTNBD030/FIP-XI PLAN). P.V. thanks DST, India for Major Research Project (SB/S1/PC-52/2012) and CSIR, India for the award of Emeritus Scientistship (Ref. No. 21(0936)/12/EMR-II).

References

1. P. M. Vignais and B. Billoud, Chem. Rev., 2007, 107, 4206–4272.
2. P. E. Siegbahn and M. R. Blomberg, Annu. Rev. Phys. Chem., 1999, 50, 221–249.
3. M. Pavlov, M. R. Blomberg and P. E. Siegbahn, Int. J. Quantum Chem., 1999, 73, 197–207.
4. L. De Gioia, P. Fantucci, B. Guigliarelli and P. Bertrand, Inorg. Chem., 1999, 38, 2658–2662.
5. M. Pavlov, P. E. Siegbahn, M. R. Blomberg and R. H. Crabtree, J. Am. Chem. Soc., 1998, 120, 548–555.
6. S. A. Vedha, R. V. Solomon and P. Venuvanalingam, Phys. Chem. Chem. Phys., 2015, 16, 10698–10707.
7. N. Argaman and G. Makov, Am. J. Phys., 2000, 68, 69.
8. E. J. Baerends and O. V. Gritsenko, J. Phys. Chem. A, 1997, 101, 5383–5403.
9. F. Neese, Coord. Chem. Rev., 2009, 253, 526–563.
10. P. Amara, A. Volbeda, J. C. Fontecilla-Camps and M. J. Field, J. Am. Chem. Soc., 1999, 121, 4468–4477.
11. P. E. Siegbahn and F. Himo, J. Biol. Inorg. Chem., 2009, 14, 643–651.
12. P. E. Siegbahn, C. R. Chim., 2007, 10, 766–774.
13. S. O. N. Lill and P. E. Siegbahn, Biochemistry, 2009, 48, 1056–1066.
14. P. E. Siegbahn, J. W. Tye and M. B. Hall, Chem. Rev., 2007, 107, 4414–4435.
15. P. Jayapal, M. Sundararajan, I. H. Hillier and N. A. Burton, Phys. Chem. Chem. Phys., 2008, 10, 4249–4257.
16. P. Söderhjelm and U. Ryde, J. Mol. Struct.: THEOCHEM, 2006, 770, 199–219.
17. S. Sumner, P. R. SöDerhjelm and U. Ryde, J. Chem. Theory Comput., 2013, 9, 4205–4214.
18. M. Kampa, W. Lubitz, M. Van Gastel and F. Neese, J. Biol. Inorg. Chem., 2012, 17, 1269–1281.
19. L. Hu, P. Soderhjelm and U. Ryde, J. Chem. Theory Comput., 2013, 9, 640–649.
20. M. Katouda and T. Nakajima, J. Chem. Theory Comput., 2013, 9, 5373–5380.
21. L. Hu, J. Eliasson, J. Heimdal and U. Ryde, J. Phys. Chem. A, 2009, 113, 11793–11800.
22. A. J. Mulholland, J. R. Soc., Interface, 2008, 5, 169–172.
23. G. Fauque, H. Peck Jr, J. Moura, B. Huynh, Y. Berlier, D. Dervartanian, M. Teixeira, A. Przybyla, P. Lespinat and I. Moura, FEMS Microbiol. Lett., 1988, 54, 299–344.
24. C. S. Baltazar, M. C. Marques, C. M. Soares, A. M. Delacey, I. A. Pereira and P. M. Matias, Eur. J. Inorg. Chem., 2011, 2011, 948–962.
25. M. C. Marques, R. Coelho, A. L. De Lacey, I. A. Pereira and P. M. Matias, J. Mol. Biol., 2010, 396, 893–907.
26. O. Gutiérrez-Sanz, M. C. Marques, C. S. Baltazar, V. M. Fernández, C. M. Soares, I. A. Pereira and A. L. De Lacey, J. Biol. Inorg. Chem., 2013, 1–9.
27. S. A. Vedha, G. Velmurugan, R. Jagadeesan and P. Venuvanalingam, Phys. Chem. Chem. Phys., 2015, 17, 20677–20686.
28. C. S. Baltazar, V. H. Teixeira and C. M. Soares, J. Biol. Inorg. Chem., 2012, 17, 543–555.
29. A. D. Becke, Phys. Rev. A, 1988, 38, 3098–3100.
30. J. P. Perdew, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 8822–8824.
31. M. Stein and W. Lubitz, Phys. Chem. Chem. Phys., 2001, 3, 5115–5120.
32. H. S. Shafaat, K. Weber, T. Petrenko, F. Neese and W. Lubitz, Inorg. Chem., 2012, 51, 11787–11797.
33. Y. Zhao and D. G. Truhlar, Acc. Chem. Res., 2008, 41, 157–167.
34. S. Grimme, A. Hansen, J. G. Brandenburg and C. Bannwarth, Chem. Rev., 2016, 116, 5105–5154.
35. C. Greco, M. Bruschi, L. De Gioia and U. Ryde, Inorg. Chem., 2007, 46, 5911–5921.
36. P. E. Siegbahn, J. W. Tye and M. B. Hall, Chem. Rev., 2007, 107, 4414–4435.
37. C. F. Macrae, P. R. Edgington, P. Mccabe, E. Pidcock, G. P. Shields, R. Taylor, M. Towler and J. V. D. Streek, J. Appl. Crystallogr., 2006, 39, 453–457.
38. K. Müller-Dethlefs and P. Hobza, Chem. Rev., 2000, 100, 143–168.
39. M. Frisch, G. Trucks, H. Schlegel, G. Scuseria, M. Robb, J. Cheeseman, G. Scalmani, V. Barone, B. Mennucci and G. Petersson, Rev. B. 01, Gaussian Inc., Wallingford CT, USA, 2009.
40. J. C. Gordon and G. J. Kubas, Organometallics, 2010, 29, 4682–4701.
41. V. Jonas and W. Thiel, J. Chem. Phys., 1995, 102, 8474–8484.
42. V. Jonas and W. Thiel, Organometallics, 1998, 17, 353–360.
43. T. Steiner, Angew. Chem., Int. Ed., 2002, 41, 48–76.
44. T. Steiner, Angew. Chem., Int. Ed., 2002, 41, 48.
45. R. E. Rosenfield Jr, R. Parthasarathy and J. Dunitz, J. Am. Chem. Soc., 1977, 99, 4860–4862.
46. J. S. Murray, P. Lane, T. Clark and P. Politzer, J. Mol. Model., 2007, 13, 1033–1038.
47. A. Upendra and S. Steve, Chem. Phys. Lett., 2011, 514, 36–39.
48. J. S. Murray, P. Lane and P. Politzer, Int. J. Quantum Chem., 2008, 108, 2770–2781.
49. S. Scheiner, J. Chem. Phys., 2011, 134, 164313.
50. S. Scheiner, Phys. Chem. Chem. Phys., 2011, 13, 13860–13872.
51. S. Scheiner, J. Chem. Phys., 2011, 134, 094315.
52. J. Contreras-García, E. R. Johnson, S. Keinan, R. Chaudret, J.-P. Piquemal, D. N. Beratan and W. Yang, J. Chem. Theory Comput., 2011, 7, 625–632.
53. S. J. Grabowski, Chem. Rev., 2011, 111, 2597–2625.

---

Footnote

† Electronic supplementary information (ESI) available: The topological properties at the BCPs of H-bonds formed between hydrophobic residues and the CO. See DOI: 10.1039/c6ra11295a

---