Metal–ligand bonding and noncovalent interactions of mutated myoglobin proteins: a quantum mechanical study

Local vibrational mode analysis

The local vibrational mode theory is briefly introduced here. Readers are encouraged to consult previous review articles for a detailed discussion of the mathematical formalism and applications. ^{22 , 23} Normal vibrational modes provide valuable electronic structure information. ^{24 , 25} However, as Wilson already pointed out in 1941, ²⁶ these modes are generally delocalized in polyatomic molecules. Consequently, the common practice of using associated normal mode stretching force constants as bond strength descriptors is questionable. In contrast, local vibrational modes, derived from normal modes, offer local vibrational mode force constants that reflect the intrinsic strength of a chemical bond or weak chemical interaction, such as the metal–ligand interactions, H-bonds, and NCIs investigated in this work.

A local vibrational mode vector a _n is defined as:

where K corresponds to the diagonal normal mode force constant matrix expressed in normal mode coordinates Q and d _n corresponds to the nth normal mode vector in internal coordinates q. Both, quantities can be found in the output of a standard normal mode calculation, offered in most modern quantum chemistry packages, ²⁷ i.e., the local mode analysis (LMA) can be applied with minimal additional computational cost. For each local mode a _n , the associated local force constants k n a describing the local vibration of a given fragment (e.g., bond, bond angle, dihedral angle, etc.) is defined as:

Additional local mode properties, such as local vibrational frequencies and intensities can be derived as well. ^{22 , 23}

We have successfully applied local vibrational mode theory to assess the strength of covalent and non-covalent interaction across the periodic table ^{22 , 23} and to characterize π–π interactions, metal–π interactions, and π–H interactions across various chemical environments utilizing specific local mode force constants, designed for this purpose. Recent applications include uranium sandwich complexes, ²⁸ metallocenes, ²⁹ and enzymes. ^{30 , 31 , 32 , 33 , 34 , 35} Additional valuable features of local mode theory include the Characterization of Normal Modes (CNM) procedure. ^{36 , 37} CNM decomposes each normal mode into local mode contributions, enhancing the interpretation of vibrational spectra for diverse systems. These range from biological structures, such as DNA base pairs, ³⁸ the exploration of the Jahn–Teller effect in metal–ligand bonding, ³⁹ to the evaluation of potential vibrational Stark-effect probes. ³⁶

In this study, we focused on local mode force constants (k ^a ) assessing metal–ligand bonding, hydrogen bonding within the distal pocket of the Schiff-base Mb, and NCIs within the active pocket of Mb. Utilizing the generalized Badger’s rule, ^{40 , 41} we transformed the local force constants into relative bond strength orders (BSO) via a power relationship BSO = A(k ^a ) ^B . Two reference molecules with known local force constants and BSO are used to obtain parameters for A and B, with the constraint that a zero value for the local force constant k ^a yields a zero BSO value. The reference molecules used in our study are presented in Table 1, where for metal–ligand bonds we designated scaled Mayer’s bond orders ^{42 , 43 , 44} instead of using the classical description of BSO = 1 for single bonds and BSO = 2 for double bonds.

We complemented LMA with Bader’s quantum theory of atoms in molecules (QTAIM) ^{45 , 46 , 47} to analyze the electron density ρ(r). Specifically, we utilized the Cremer-Kraka criterion ^{48 , 49} to characterize the covalent and electrostatic nature of bonds based on the local energy density H(r) = G(r) + V(r), where the kinetic density is G(r) and the potential energy is V(r). According to this criterion, the bond character between two atoms (A and B) at the bond critical point can be determined as follows: a negative H(r _b ) indicates predominantly covalent character, while a positive H(r _b ) signifies predominantly electrostatic character. In addition, we calculated the charge difference of the metal and oxygen utilizing natural population analysis (NPA) to further investigate the M–L bond interactions. ⁵⁰

Computational methods

We first began by taking PDB structures (1V9Q for Mn(III) salophen, ¹⁹ 1UFJ for Fe(III) salophen, ¹⁸ and 1J3F for Cr(III) salophen ¹⁹ ) to model salophen. PDBs with salophen were taken as the starting point and edited manually to model salen structures. Missing amino acids within the protein chain were added with Modeller in ChimeraX. ^{51 , 52} The sequences of missing amino acids can be found in the Supporting Information, Section 1. A similar procedure was conducted for cyanide ligands, using oxygen from water coordination as the starting point. For systems with mutation of H64D, the mutation was applied and the best score was selected with Dynameomics rotamer in ChimeraX. ^{52 , 53} All prosthetic groups with metals and water/cyanide coordination were modeled with Metal Center Parameter Builder (MCPB.py) to prepare the protein structure. ⁵⁴ This approach ensures that molecular modeling of the metal in the protein is accounted for appropriately. Starting from pre-determined models that include the metal, Schiff base, and proximal/distal histidine (H93, H64) up to the β-carbons, QM optimization, and frequency calculations were done with ωB97X-D functional ⁵⁵ and Pople’s 6-31G(d,p) basis set ^{56 , 57} as implemented with Gaussian16. ⁵⁸ Bonded parameters were obtained from Seminario method. ⁵⁹ For Cr, Mn, and Fe ions, as well as coordinating ligands, the Merz–Singh–Kollman atomic partial charges ⁶⁰ were obtained with the aforementioned model chemistry and the restrained electrostatic potential (RESP) approach. ⁶¹ The coordinating ligands’ backbone (heavy and hydrogen) atoms were held to standard force field values. Metal center parameters were developed for use with the ff19SB AMBER forcefield. ⁶² After modeling the metal center, hydrogens were added with H++ server, ^{63 , 64 , 65} and the proteins were neutralized with Cl-atoms, then a water sphere of OPC waters was added 16 Å from the metal center utilizing tleap from AMBER. ⁶⁶ Two protonation forms were modeled for the distal histidine: epsilon form (H64ϵ) and delta form (H64δ). After appropriate modeling of the protein’s metal center and corresponding waters in AMBER, H’s were minimized for 2000 steps using the steepest descent. TAOpackage was utilized to prepare minimized protein structures as inputs for QM/MM optimization. ⁶⁷ The QM region incorporates the salophen/salen structure coordinated to the water and the proximal (H93) and distal (H64) histidines (average ∼70 atoms), while the remaining portion of the protein is considered the MM region (average ∼3204 atoms). All proteins were first optimized with mechanical embedding employing the ωB97X-D/6-31G(d,p)/AMBER level of theory, followed by electronic embedding with scaled charges (scalecharge = 333 300) to account for the polarization of the MM region in the QM wavefunction. ^{55 , 56 , 57} For both optimizations, the self-consistent field (SCF) was set to quadratic convergence, and optimized with quadmacro. ^{68 , 69} To account for potential spin state effects, we evaluated both high-spin and low-spin states for representative systems. These tests were conducted in the gas phase and with a few QM/MM systems, and the energetically preferred spin states for the water ligand were used for further analysis. For the CN⁻ ligand, depending on the metal, it can either evoke a high-spin or low-spin state; therefore, for consistency, we have calculated all systems with CN⁻ in the high-spin state. Detailed spin-state energetics and structural comparisons are provided in the Supporting Information (Section 2). Frequency calculations ensured no imaginary frequencies to confirm the protein structures were at a local minima. LMA was done with the LMODEA23 package to study metal–ligand bonding, hydrogen bonding, and NCIs. ⁷⁰ QTAIM analysis was conducted with AIMALL, ⁷¹ and NBO analysis was done with NBO version 7.0. ⁷²

Metal–ligand bonding

This section explores the trends and nature of M–L bonding in Cr, Mn, and Fe systems, using local force constants, QTAIM energy densities, and natural population analysis (NPA) charge (as summarized in Table 2) to elucidate the relationship between structure and bonding characteristics. We first analyze trends in metal–ligand bonding for Cr, Mn, and Fe systems, focusing on the influence of ligand type (water, OH₂, and cyanide, CN⁻, herein labeled CN), prosthetic group type (salphen, S1 and salen, S2) as well as distal histidine (H64ϵ, H64δ, and H64D). Figures 2(a)–(d), 3(a)–(d), and 4(a)–(d) display the different trends of the metal–ligand bond(s) concerning the bond strength (k ^a , mdyn/Å).

Fig. 2:

Properties of Cr–L bonds (L = OH₂, CN) in Cr(III) Mb Schiff bases and corresponding gas phase models. (a) Relation between Cr–L bond length (Å) and local force constant k ^a (Cr–L). R ² values are calculated based on power regression fit. Spearman correlation coefficient ρ is also displayed. (b) BSO n (Cr–L) calculated from local force constant via generalized Badger rule. (c) Relation between energy density H(r) and local force constant k ^a (Cr–L). (d) Relation between atomic charge difference Δe and local force constant k ^a (Cr–L).

Fig. 3:

Properties of Mn–L bonds in Mn(III) Mb Schiff bases and corresponding gas phase models. (a) Relation between Mn–L bond length (Å) and local force constant k ^a (Mn–L). R ² values are calculated based on power regression fit. Spearman correlation coefficient ρ is also displayed. (b) BSO n (Mn–L) calculated from local force constant via generalized Badger rule. (c) Relation between energy density H(r) and local force constant k ^a (Mn–L). (d) Relation between atomic charge difference Δe and local force constant k ^a (Mn–L).

Fig. 4:

Properties of Fe–L bonds in Fe(III)Mb Schiff bases and corresponding gas phase models. (a) Relation between Fe–L bond length (Å) and local force constant. R ² values are calculated based on power regression fit. Spearman correlation coefficient ρ is also displayed. (b) BSO n (Fe–L) calculated from local force constant via generalized Badger rule. (c) Relation between energy density H(r) and local force constant k ^a (Fe–L). (d) Relation between atomic charge difference Δe and local force constant k ^a (Fe–L).

For the chromium systems, Fig. 2(a) illustrates the correlation between bond length and strength, while Fig. 2(b) displays the local force constant as a function of BSO. Despite the correlation of bond length via power regression (R ² = 0.7041) and local force constant in Fig. 2(a) appearing to have a weak correlation, we have also shown the Spearman correlation coefficient, ρ, to be consistent with the generalized Badger rule, which states that a shorter bond equates to a stronger bond strength, ⁴¹ where the negative sign of ρ reflects the expected inverse nature of the relationship. These analyses reveal that the strongest bonds are observed in systems where the distal histidine is mutated to aspartate (H64D), with the salophen (S1) (Cr-S1-H64D, 2.173 mdyn/Å) exhibiting slightly stronger bond compared to the salen (S2) (Cr-S2-H64D, 1.836 mdyn/Å), despite being similar bond lengths (1.917 Å and 1.940 Å, respectfully). Another noteworthy trend is the effect of the type of histidine on Cr–OH₂ bonding. In both S1 and S2, the epsilon form of histidine (Cr-S1-H64 ϵ and Cr-S2-H64 ϵ) weakens the Cr–OH₂ bond (0.631 mdyn/Å and 0.664 mdyn/Å) relative to the gas-phase structure (Cr-S1-Gas, 1.113 mdyn/Å and Cr-S2-Gas, 1.125 mdyn/Å), while the delta form (Cr-S1-H64 δ, 1.261 mdyn/Å and Cr-S2-H64 δ, 1.268 mdyn/Å) strengthens it. A similar pattern is observed for Cr–CN bonds: gas-phase structures exhibit the strongest bonds (Cr-S1-CN-Gas, 1.881 mdyn/Å and Cr-S2-CN-Gas, 1.851 mdyn/Å), with the S1 slightly stronger than S2 (average of S1 with CN ligand is 1.829 mdyn/Å compared to S2 average 1.795 mdyn/Å). For both Cr–OH₂ and Cr–CN bonds, the H64δ form results in stronger bonds than the H64ϵ form, highlighting the influence of prosthetic group and histidine type on Cr–L bonding. Furthermore, Fig. 2(b) confirms that all Cr–OH₂ and Cr–CN bonds have bond strength orders of less than one, indicating that these interactions are characterized as less than a formal single bond.

We now examine the nature of these bonds by analyzing energy density and natural population charges. Figure 2(c) compares the energy density (H(r), Hartree/ų) with the local force constant for the Cr–L bonds. According to the Cremer–Kraka criterion, a negative energy density indicates covalent bonding, while a positive or near-zero energy density suggests electrostatic (noncovalent) interactions. ^{73 , 74} This analysis reveals that Cr–OH₂ gas-phase structures, as well as the H64ϵ and H64δ forms of both S1 and S2, are characterized as electrostatic Cr–OH₂ bonds. This is reflected in the weaker local force constants, which show a slight correlation between bond strength and energy density. However, the correlation coefficient values are relatively low, indicating a weak correlation. The outliers to this trend are the Cr–OH₂ bonds in the H64D mutation, where proton transfer from water to aspartate induces a shift toward covalent bonding. This change is reflected in stronger local force constants and negative energy densities, signaling a more stabilized bond. For Cr–CN bonds, similar trends are observed, where gas-phase structures exhibit stronger bonds with greater covalency, while the H64ϵ form in both S1 and S2 results in weaker Cr–CN bonds and less covalent character. Cyanide ligands are known to be stronger ligands in the spectrochemical series, primarily due to their π-backbonding capabilities, which enhance their covalent interactions with metal centers. This stronger interaction is consistent with the observed trend that cyanide ligands form more covalent bonds with chromium compared to water ligands. However, as with the Cr–OH₂ systems, the correlation between bond strength and energy density for the Cr–CN bonds remains weak, as reflected in the low correlation coefficients. Further insight is provided by the charge difference between chromium and the ligands, as shown in Fig. 2(d). The Cr–OH₂ bond exhibits a charge difference ranging from 2.781 to 2.927, indicating a certain degree of charge transfer between chromium and the water ligand that is modulated by the electronic structure, based on the prosthetic group and histidine type. In contrast, the charge difference for the Cr–CN bond is significantly smaller (below 2.0), reflecting a more balanced charge distribution, which is characteristic of stronger metal–ligand interactions, underscoring the nature of cyanide as a stronger, more covalently bound ligand compared to water. It is important to note that NBO charge differences do not provide a strong predictor of bond strength, suggesting that they should be considered as one factor among others.

Figure 3(a) and (b) show the relationship between bond length and bond strength for Mn–L interactions, as well as the BSO plot with local force constants, respectively. There is a higher correlation (R ² = 0.9744) of bond length and local force constant, (coupled with a strong negative ρ coefficient), with the general trends being similar to those in the chromium system. Specifically, the H64ϵ form results in a weaker and longer Mn–OH₂ bond (Mn-S1-H64 ϵ, 0.127 mdyn/Å, 2.623 Å and Mn-S2-H64 ϵ, 0.219 mdyn/Å, 2.525 Å) compared to the gas-phase structures (Mn-S1-Gas, 0.300 mdyn/Å, 2.405 Å and Mn-S2-Gas, 1.023 mdyn/Å, 2.104 Å), while the H64δ form strengthens the Mn–OH₂ bond for the S1 (Mn-S1-H64 δ, 0.660 mdyn/Å). However, the S2 in the H64δ form (Mn-S2-H64 δ, 0.407 mdyn/Å) exhibits a weaker Mn–OH₂ bond than its gas-phase counterpart. This may be due to the lengthening of the Mn–OH₂ bond as it goes from gas phase 2.104 Å to 2.290 Å in the S2 of H64δ. In contrast to the chromium system, where the strongest metal–ligand interactions are observed for the H64D mutation with Cr–OH₂ bonds, the strongest Mn–L bonds are found with cyanide ligands. For manganese, the strongest bond occurs in the gas phase with the S1 (Mn-S1-CN-Gas), followed by the H64δ form of S1 (Mn-S1-CN-H64 δ) and the H64ϵ form of S2 (Mn-S2-CN-H64 ϵ). Figure 3(b) follows the trends observed in chromium, with the BSO of Mn–L bonds consistently showing values characteristic of interactions that are less than a single bond.

Unlike the chromium systems, Fig. 3(c) shows that there is less correlation between the water and cyanide ligands in terms of the local force constant. All of the Mn–L bonds exhibit covalent character, with stronger bonds demonstrating greater covalency than weaker bonds. In Fig. 3(d), the NBO charge difference (Δe) reveals a correlation between charge distribution and the local force constant, although there are two notable outliers for the Mn–OH₂ bond (Mn-S2-H64D and Mn-S2-Gas). As observed previously, the water ligands show a greater charge difference between manganese and the oxygen of the water, while for cyanide ligands, the charge difference decreases significantly. The charge difference for cyanide ligands typically ranges from 1.992 to 2.095, suggesting a more balanced charge distribution compared to water ligands.

Lastly, for iron, Fig. 4(a) displays the relationship between bond length and bond strength for Fe–L interactions, while Fig. 4(b) shows the BSO plot with the local force constant for Fe–L bonds. Here, the correlation remains for bond length and local force constant (R ² = 0.9175), with the observed trends being similar to those seen for chromium and manganese. Specifically, the distal histidine (H64ϵ) form results in the weakest and longest bonds for both salophen (S1) and salen (S2) (represented by right triangles, Fe-S1-H64 ϵ, 0.030 mdyn/Å, 2.683 Å and Fe-S2-H64 ϵ, 0.012 mdyn/Å, 2.677 Å). We have previously investigated aquamet-Mb, a water-ligand heme Mb protein, where it was determined that H64δ form was a stronger Fe–OH₂ bond compared to H64ϵ form. ^{30 , 37} This shows that the mutation of distal histidine to aspartate (H64D) as well as the inclusion of cyanide ligands, leads to stronger Fe–L bonds in both protonation states of the distal histidine. The strongest Fe–L bond is observed in the H64D mutation with OH₂ ligand, with the salophen (Fe-S1-H64D, 1.848 mdyn/Å) showing slightly stronger Fe–OH₂ bond strength than the salen (Fe-S2-H64D, 1.621 mdyn/Å). Figure 4(b) further supports these trends, with Fe–L bonds consistently characterized by bonds less than a single bond, as seen in chromium and manganese systems.

As shown in Fig. 4(c), the Fe–OH₂ bonds in the H64ϵ forms for both S1 and S2 (Fe-S1-H64 ϵ and Fe-S2-H64 ϵ) exhibit almost exclusively electrostatic character, with values indicating minimal covalent interaction (−0.0087 Hartree/ų and -0.0089 Hartree/ų, respectively). In contrast, the remaining structures demonstrate greater covalent character. This trend is consistent with Fe–OH₂ ligands showing less covalency compared to Fe–CN bonds. We also see a trend of reduced covalency that parallels the trend in Lewis acidity (Fe³⁺ > Mn³⁺ > Cr³⁺), where higher Lewis acidity enhances the metal’s ability to form stronger covalent bonds. As the metal’s Lewis acidity decreases, the bonding becomes more electrostatic, as with Cr–OH₂ bonds, and can be reflected in weaker interactions with water ligands. This observation aligns with the behavior of the distal histidine residue: in the H64ϵ tautomer, the imidazole group acts as a Lewis acid, withdrawing electron density from the water ligand and thereby reducing its Lewis basicity, and generally leads to weaker, more electrostatic metal–ligand interactions. Conversely, the H64δ tautomer acts as a Lewis base, enhancing the water ligand’s donor ability and favoring stronger, more covalent bonding with the metal center. When analyzing the NBO charge difference (Δe) in Fig. 4(d), we observe that the charge difference for Fe–OH₂ bonds is greater than that for the corresponding Cr and Mn systems, suggesting a stronger electrostatic component in the Fe–OH₂ interactions. The Fe–CN ligands, however, show less variation in both force constant and charge difference, indicating more consistent bonding characteristics across different prosthetic groups and protonation states.

Noncovalent interactions

Due to the stability of these Schiff Base prosthetic groups in apomyoglobin protein via experimental characterization, we sought to understand the interplay of NCIs that transpire at the atomistic level. We explored two avenues of NCIs that occur: hydrogen bonding at the distal site (H64 or H64D), and π–π or π–H interactions that take place at the second coordination sphere and active site of myoglobin.

Hydrogen bonding

We first investigated the role of hydrogen bonding (H-bond) between the distal histidine (in both forms of protonation state) or its mutation to aspartate (H64D), and the ligands (both water and cyanide). Table 3 exhibits systems with H-bond’s distance (d, Å), local force constant (k ^a , mdyn/Å), BSO (H…N/O), and energy density (H(r), Hartree/ų). The correlation between bond length and local force constant is consistent with the Badger rule, where a stronger bond indicates a shorter bond (Fig. 5(a), (c) and (e)). Several types of H-bonds are observed depending on the distal site of myoglobin. When the histidine is in the H64ϵ form, it acts as the H-bond donor, while either the oxygen from water or the nitrogen from cyanide serves as the H-bond acceptor. In the H64δ form, water acts as the H-bond donor, with nitrogen from cyanide accepting the hydrogen. However, no H-bond donor is present for cyanide ligands. In this case, the histidine undergoes a distortion that results in a very weak interaction between a proton from carbon in histidine and the nitrogen in cyanide (see Figs. S1–S5 in Supporting Information for a depiction of interactions). Finally, for the H64D mutation, proton transfer occurs from one of the hydrogens of water to the oxygen of aspartate, creating an H-bond with aspartate as the donor and water as the acceptor. The H-bond distances range from 1.4 to 2.9 Å, with an average distance of 1.96 Å.

Table 3:

Table showing the distance, d between the H…O/N bond (in Å), local force constant, k ^a (mdyn/Å), bond strength order, BSO (H…O/N), and energy density H(r) (Hartree/ų).

Name	d	k a	BSO (H…O/N)	H(r)
Cr-S1-H64ϵ	2.026	0.188	0.272	−0.006
Cr-S2-H64ϵ	1.994	0.177	0.267	−0.005
Cr-S1-H64δ	1.645	0.317	0.323	−0.063
Cr-S2-H64δ	1.651	0.328	0.327	−0.055
Cr-S1-H64D	1.553	0.304	0.318	−0.064
Cr-S2-H64D	1.477	0.228	0.290	−0.151
Cr-S1-CN-H64ϵ	2.278	0.092	0.215	0.010
Cr-S2-CN-H64ϵ	2.275	0.074	0.200	0.009
Cr-S1-CN-H64δ	2.544	0.055	0.182	0.005
Cr-S2-CN-H64δ	2.587	0.057	0.184	0.006
Mn-S1-H64ϵ	1.948	0.215	0.284	−0.006
Mn-S2-H64ϵ	2.035	0.157	0.256	−0.006
Mn-S1-H64δ	1.651	0.394	0.347	−0.054
Mn-S2-H64δ	1.769	0.307	0.320	−0.013
Mn-S1-H64D	1.448	0.223	0.288	−0.203
Mn-S2-H64D	1.455	0.377	0.342	−0.167
Mn-S1-CN-H64ϵ	2.031	0.173	0.265	0.005
Mn-S2-CN-H64ϵ	2.045	0.123	0.237	0.005
Mn-S1-CN-H64δ	2.584	0.052	0.178	0.006
Mn-S2-CN-H64δ	2.485	0.100	0.221	0.008
Fe-S1-H64ϵ	1.902	0.229	0.290	−0.004
Fe-S2-H64ϵ	1.904	0.180	0.268	−0.030
Fe-S1-H64δ	1.713	0.306	0.319	−0.004
Fe-S2-H64δ	1.708	0.307	0.320	−0.030
Fe-S1-H64D	1.617	0.251	0.299	−0.031
Fe-S2-H64D	1.521	0.238	0.294	−0.096
Fe-S1-CN-H64ϵ	2.206	0.085	0.210	0.008
Fe-S2-CN-H64ϵ	2.083	0.129	0.240	0.005
Fe-S1-CN-H64δ	2.953	0.023	0.137	0.008
Fe-S2-CN-H64δ	2.440	0.073	0.199	0.005

Fig. 5:

Properties of hydrgoen bonding in M(III) Schiff bases. (a) Cr hydrogen bonding correlation between bond length and local force constant. R ² values are calculated based on power regression fit. Spearman correlation coefficient ρ is also displayed. (b) Cr hydrogen bonding BSO (c) Mn hydrogen bonding correlation between bond length and local force constant. (d) Mn hydrogen bonding BSO. (e) Fe hydrogen bonding correlation between bond length and local force constant. (f) Fe hydrogen bonding BSO ωB97X-D/6-31G(d,p)/AMBER level of theory.

Figure 5(b) shows the bond strength order (BSO) (H…N/O) of the chromium systems in relation to their local force constant. The strongest H-bond interactions in the chromium systems occur with Cr–OH₂ ligands. Specifically, Cr-S2-H64 δ (0.328  mdyn/Å and 1.651 Å) and Cr-S1-H64 δ (0.317 mdyn/Å and 1.645 Å) exhibit the strongest H-bond interactions. Comparatively, the H64ϵ form of histidine creates less favorable H-bonding geometries, leading to longer bonds and weaker interactions, as seen in Cr-S1-H64 ϵ (0.188  mdyn/Å and 2.206 Å) and Cr-S2-H64 ϵ (0.177  mdyn/Å and 2.026 Å). The mutation of the distal histidine to aspartate (H64D) creates a moderate H-bond interaction. Additionally, the presence of cyanide (Cr–CN) ligands weakens the H-bond interactions, with the weakest observed in Cr-S1-CN-H64 δ (0.055 mdyn/Å and 2.544 Å) and Cr-S2-CN-H64 δ (0.057 mdyn/Å and 2.587 Å).

For manganese, the trends differ slightly, as shown in Fig. 5(d). The strongest H-bonding interactions occur in Mn-S1-H64 δ and Mn-S2-H64D, demonstrating the adaptability of how both the metal–ligand bond and secondary interactions, such as H-bonding, influence the overall interaction, as captured by the local force constants. Mn-S2-H64 δ and Mn-S1-H64D exhibit moderate strength, followed by Mn-S1-H64 ϵ and Mn-S1-CN-H64 ϵ. These results suggest that, for water ligands, there is a slight preference for H64D and H64δ systems, though the trend is dependent on the prosthetic group type: preference for S1 when histidine is present and preference for S2 when H64D is present. The weakest H-bond occurs with the cyanide ligands, in Mn-S1-CN-H64 δ (0.052  mdyn/Å and 2.584 Å).

Following the manganese systems, the iron systems shown in Fig. 5(f) exhibit the strongest H-bond interactions in Fe-S2-H64 δ (0.307  mdyn/Å and 1.708 Å) and Fe-S1-H64 δ (0.306  mdyn/Å and 1.713 Å), followed by Fe-S1-H64D. Moderate H-bond strength is observed in Fe-S2-H64D and Fe-S1-H64 ϵ. The weakest H-bond occurs in Fe-S1-CN-H64 δ (0.023 mdyn/Å and 2.953 Å). Overall, the trends show that CN ligands exhibit weaker H-bonds overall, with H64δ being the weakest. Unlike the water ligand environment, where hydrogen from the water acts as the H-bond donor when switching from H64ϵ to H64δ, the CN ligand does not have an H-bond donor, and as discussed previously, the distal histidine distorts to have hydrogen from a nearby carbon form a weak nonclassical hydrogen bond (seen in Fig. S4). This can also be realized through energy density H(r) in Table 3, where all CN H-bonds are positive, indicating electrostatic character, unlike the rest of the H-bonds having negative covalency character.

H-bonding visualized by Characterization of Normal Modes

Since our investigation focused on the elucidation of differences in H-bonding at the distal histidine of the mutated Mb, we utilized in addition to the comparison of H-bond strengths CNM as an important tool to uncover sensitive differences in the H-bond normal mode composition. It is important to note that the CNM procedure requires a complete non-redundant and chemically meaningful set of N _vib local mode parameters, with N _vib = (3N − 6) for nonlinear molecule and (3N − 5) for a linear molecule being composed of N atoms. For QM/MM systems where the total number of QM and MM atoms can easily reach several thousand, this can be daunting task. To address this challenge, our group has developed two tools to streamline the CNM procedure for large systems. First, the Generalized Subsystem Vibrational Analysis (GSVA) projects the full N _vib vibrational space of a large system, e.g., a QM/MM system into a reduced, user-defined subspace, such as the QM region of the QM/MM system. ⁷⁵ Next, LModeAGEN (integrated in LModeA) which is based on graph theory, automatically generates a complete set of non-redundant, chemically meaningful local mode parameters. ^{76 , 77}

As an example, CNM was implemented for Cr-S1-H64 ϵ and Cr-S1-H64 δ, both containing 3264 atoms, resulting in a total vibrational mode count (N _vib ) of 9786. By restricting the analysis to the 72 QM atoms, N _vib is reduced to 204 vibrational modes. Figure 6(a) and (b) illustrate a variety of normal vibrational modes, highlighting significant contributions from either the Cr–OH₂ or Cr–N(H93) interactions, typically observed at lower frequencies, as well as stretching modes associated with the H-bonding distal histidine region. Figure 6(a), which depicts the CNM for Cr-S1-H64 ϵ, shows a strongly localized normal mode at 3663 cm⁻¹ for the distal histidine N–H stretch, with 96 % local mode character. At the lower frequency of 61 cm⁻¹, significant delocalization is observed, with only 16 % local mode character attributed to the H-bond (O…H). The O–H asymmetric stretches for water, occurring at 3758 cm⁻¹ and 3870 cm⁻¹, display varying degrees of delocalization, with local mode contributions of approximately 62 % and 36 %, respectively. In contrast, Fig. 6(b), which represents the CNM for Cr-S1-H64 δ, highlights the distinct nature of H-bonding associated with different protonation states of histidine. One O–H bond participating in an H-bond interaction with histidine results in a completely localized O–H stretch at 3921 cm⁻¹, with 99 % local mode character. Additionally, at 2702 cm⁻¹, the O–H bond involved in H-bonding to the distal histidine exhibits 92 % local mode character, with 7 % of the local mode contribution arising from the H-bond interaction (N…H). Utilizing CNM as a powerful tool, these findings highlight the pivotal role of protonation states and H-bonding interactions in shaping the vibrational properties and mode localization within the protein system.

Fig. 6:

Decomposition of selected normal vibrational modes into stretching, bending, and torsional local modes involving H-bond between water and distal histidine in myoglobin (a) in Cr-S1-H64 ϵ and (b) in Cr-S1-H64 δ.

Noncovalent local force constants

NCIs and their role in modifying metalloproteins have been extensively discussed in previous reviews and articles. ^{21 , 78 , 79 , 80} In the active pocket, or secondary coordination sphere, a variety of weak interactions occur with the heme prosthetic group in typical myoglobin proteins. ^{79 , 81} We have previously investigated π–π interactions for mutated tryptophan Mb’s utilizing a unique local force constant that takes the geometric center point of a ring and projects the stretching component utilizing Cremer-Pople’s curvilinear coordinates. ³⁰ Additionally, studies validating the characterization of NCIs within a QM/MM framework demonstrated that the AMBER force field captures dispersion interactions with lower mean average error than other commonly used force fields. ⁸² However, these findings were based on cumulative energetic properties. We propose that our local mode approach offers a more direct and computationally efficient means of capturing these weak interactions, without sacrificing accuracy. This method could enable the efficient study of secondary interactions within metalloprotein systems, providing new insights into their role in structural and functional modulation.

To capture the range of weak interactions that could occur, we have narrowed the possibility of interactions to be within 4 Å of the active size, utilizing VMD software ⁸³ to quantify the amino acids of interest. We then only considered two types of interactions: π–π and π–H interactions that are interacting with the Schiff base prosthetic group. By doing this, we can capture the direct NCIs local force constant that might affect the prosthetic group’s geometry. Overall three amino acids will be discussed: F43, V68, and H97. Their relative placement concerning the active pocket is shown in Fig. 7. It is important to note that certain interactions vary depending on the orientation of the Schiff base prosthetic group’s ring system. To address this, Fig. 8 provides a schematic with labels (A–F) representing the types of NCIs studied in myoglobin (Mb). Overall, six types of NCIs were identified across the 30 proteins analyzed. A brief description of each type is as follows: H97-A is either a π–π interaction for S1 or a π–H interaction for S2, since it lacks the additional ring. V68-B and V68-F are both π–H interactions, while F43-A/E are both π–π interactions, with F43-A occurring for S1 and F43-E occurring for S2. F43-C and F43-D are both π–π interactions regardless of whether it is an S1 or S2.

Fig. 7:

Noncovalent interactions of interest that occur within Schiff base myoglobin.

Fig. 8:

Labeling of rings within Schiff base salophen.

A preliminary analysis showed no obvious correlation between bond length and bond strength, showcasing that in this scenario the Badger rule does not apply (Fig. S6 in Supporting Information). It is also noteworthy to mention that due to the nature of the projection procedure (where we calculate utilizing Cremer’s mean ring plane ⁸⁴ in the z direction ³⁰ ) to calculate the unique local force constant, it creates a shorter local bond length than if one were to use normal parameters in LMODEA. A complete table of all NCI parameters calculated is shown in Table S6 of the Supporting Information, with the main focus presented here utilizing the Sz stretching (indicated by index 23 in LMODEA). The greatest interactions usually occur with H97-A, with weaker interactions for S2 (since it is a π–H interaction rather than a π–π interaction). However, there is a lot of variability, showcasing that every protein environment is different, which is captured by local force constants. While local mode force constants derived from harmonic analysis provide insight into the intrinsic vibrational stiffness of noncovalent interactions, it is important to emphasize the limitations of harmonic potentials in accurately capturing weak dispersion interactions and cooperative effects. These complexities likely contribute to the lack of correlation between local force constants and bond distances observed here. We therefore outline the potential for future work to investigate these aspects in greater detail, possibly through comparisons with more advanced models or pairwise dispersion energy calculations.

To better understand the influence of the collective 6 NCI interactions that occur, we have collected and plotted the average NCI local force constants in Fig. 9. From here we can make observations on the collective trends that affect the interaction of NCI from the type of prosthetic group, ligand, distal histidine protonation form, as well as the metal. Looking first at chromium, the highest average NCI local force constant is Cr-S2-H64 ϵ while the lowest average is Cr-S2-H64 δ. The general trend is that the S2 and H64ϵ create stronger NCIs, while for the H64D mutation, the S2 decreases in strength compared to S1. For CN ligands, minimal changes occur within the NCIs, indicating that the effect of the ligand does not induce a change in the secondary coordination sphere as much as the change of the distal histidine for chromium systems. For manganese systems the highest average is Mn-S1-H64D, with the lowest coming from Mn-S2-CN-H64 ϵ. No general trend is observed for the prosthetic group in H64ϵ or H64δ form, however there is a significant decrease in average local force constant going from S1 to S2 for H64D mutations. For CN ligands the trend is flipped between H64ϵ and H64δ, where S2 has greater local force constant than S1 with H64δ, and S1 has slightly greater local force constant than S2 with H64ϵ. Lastly, for iron the highest average local force constant is Fe-S2-H64D, with the lowest being Fe-S1-CN-H64 ϵ. There is not much of a difference going from H64ϵ to H64δ for S1, although there is an increase of local force constant going from the same direction for S2. It appears that for H64D mutation S2 has a higher average force constant than S1. For CN ligands S2 has a higher average with H64ϵ, whereas S1 has a higher average with H64δ.

Fig. 9:

Average local force constants for noncovalent interactions. In this plot, HIE indicates H64ϵ, while HID indicates H64δ.