Exploring the noncovalent chemistry in the catalytic centers of snake venom metalloproteinases: a computational perspective

Abstract

Snakebite envenoming is a major neglected tropical disease, in which snake venom metalloproteinases (SVMPs) are key drivers of local tissue damage, microvascular disruption, and hemostatic imbalance. Closely related SVMPs can nevertheless range from highly hemorrhagic to essentially non-hemorrhagic while sharing a conserved Zn(II) catalytic site, and the structural causes of this divergence (particularly the role of the conserved methionine-containing turn (MET-turn) beneath the metal center) remain unclear. Here, we have combined a Protein Data Bank (PDB) inspection with theoretical calculations at the PBE0-D3/def2-TZVP level to investigate the noncovalent interactions stabilizing the MET residue over the Zn(II) catalytic site in both hemorrhagic and non-hemorrhagic P-I and P-III SVMPs. Additionally, the results were analyzed using several computational tools, such as the Quantum Theory of Atoms in Molecules (QTAIM) and the Non Covalent Interaction plot (NCIplot) methodologies, revealing and quantifying S lone pair–π, CH–π and Zn-associated σ-hole (Spodium bond-like) interactions. We believe the results presented herein will be useful to those scientists devoted to bioinorganic chemistry and rational drug design as well as to expand the current knowledge of SpBs among the chemical biology community.

---

Introduction

Snakebite envenoming is now recognized as a major yet long-neglected global health problem, causing tens of thousands of deaths and a large burden of disability each year, predominantly in rural communities of Africa, Asia, and Latin America.^1,2 Recent estimates suggest that 81 000–138 000 people die annually from venomous snakebites, with many more suffering permanent physical and socioeconomic sequelae.^3,4 In response, the World Health Organization (WHO) has classified snakebite envenoming as a neglected tropical disease and launched a global strategy aiming to reduce associated deaths and disabilities by 2030.^4,5 Among the multiple protein families that compose snake venoms, snake venom metalloproteinases (SVMPs) stand out as key determinants of local tissue damage, hemorrhage, and hemostatic disturbances, especially in viperid venoms.^6–8

SVMPs are part of a structurally and functionally diverse family of Zn-dependent endopeptidases that cleave extracellular matrix (ECM) components, basement membrane proteins, and a range of plasma substrates, thereby contributing to blistering, myonecrosis, coagulopathy, and particularly microvascular damage and hemorrhage.^9,10 Based on their domain architecture, SVMPs are classically classified into P-I, P-II, and P-III classes. More in detail, from an structural point of view P-I enzymes contain only the typically metalloproteinase domain; P-II toxins add a disintegrin domain; and P-III toxins include disintegrin-like and cysteine-rich domains, while forming higher-order oligomeric assemblies in certain cases.^7,11 This increasing structural complexity is related to a broadening of biological activities, since many P-III SVMPs are strongly hemorrhagic with high affinity for microvessels, whereas numerous P-I enzymes display weaker hemorrhagic effects and are more prone to fibrinogenolytic or procoagulant activities.^8,9 Therefore, understanding the structural determinants that regulate individual SVMPs hemorrhagic or non-hemorrhagic activity remains central for understanding venom pathology and for guiding the design of novel antivenoms and small-molecule inhibitors.

At the core of SVMP function lies a conserved metalloproteinase catalytic domain belonging to the metzincin superfamily of Zn-dependent endopeptidases. Metzincins are defined by a characteristic Zn-binding consensus sequence, which forms part of an extended α-helix and provides three HIS residues coordinated to the catalytic Zn(II) ion, along with a conserved “Met-turn” located directly underneath the metal site.^12,13 In this context, the analysis of high-resolution crystal structures of SVMPs and related metzincins revealed a distorted tetrahedral Zn(II) center coordinated by the N atoms of three vicinal HIS residues and a water (or hydroxide) molecule, with a nearby GLU residue acting as base/acid to activate the nucleophilic water during peptide-bond hydrolysis.¹⁴

A particular structural feature of metzincins is the MET-turn: a protein backbone 1,4-turn located immediately to the Zn(II) site containing a conserved MET residue. This MET residue inserts as a “plug” into the conserved hydrophobic pocket beneath the catalytic site, contributing to its supramolecular architecture.^15,16 In this regard, structural comparisons across metzincins showed that, despite substantial sequence differences, the position and orientation of the Met-turn residues are remarkably conserved.¹⁷ In addition, site-directed mutagenesis in matrix metalloproteinases (MMPs) and related enzymes has demonstrated that replacing this MET residue often leads to a reduced catalytic efficiency, altered thermal stability, or changes in substrate specificity, even when the global protein fold remains nearly intact.¹⁸ Altogether, these studies indicate a dual structural and functional role played by the MET-turn in (i) modulating enzymatic activity and (ii) stabilizing the Zn-binding environment. In fact, HIS⋯MET interactions have been recently proposed to control the coordination in a metal-α-helix peptide framework.¹⁹

Aside from this, this MET residue also participates in a rich network of noncovalent interactions (NCIs) that likely fine-tune the balance between structural stability and functionality in SVMPs (see Fig. 1 below). More specifically, it can be involved in (i) sulphur lone pair–π (lp–π)^20–22 contacts and (ii) CH–π interactions^23–25 involving the MET –CH₃ group and the metal coordinated HIS residues and (iii) CH⋯Zn interactions, that may be understood as Spodium Bond-like interactions (SpBs).^26,27

Fig. 1 X-ray crystal structure of 1BUD²⁸ showing the Zn(II) center alongside the MET residue. A schematic representation of the NCIs involving the MET –S–CH₃ group and the Zn(II) coordination complex are shown inside the circular part of the figure.

Given this tight balance of NCIs, subtle rearrangements in the interaction network (e.g. changes in the packing of the MET-turn against the Zn(II) site or in its contacts with the metal coordinated HIS residues) may therefore influence substrate recognition, proteolytic efficiency, and ultimately whether a given SVMP exhibits predominantly hemorrhagic or non-hemorrhagic behavior. In this regard, while molecular simulations^29–32 and structure–function studies of hemorrhagic and non-hemorrhagic SVMPs have hinted at differences in flexibility and packing around the MET-turn region,^9,33 a detailed quantum-mechanical (QM) description of the underlying noncovalent interactions remains scarce in the literature.

To tackle this point, it is advantageous to integrate structural information from the Protein Data Bank (PDB)³⁴ with QM methods capable of evaluating the subtle balance of forces that stabilize the MET-turn environment. In this work, we have combined the analysis of high-resolution crystal structures involving both hemorrhagic and non-hemorrhagic SVMPs with density functional theory (DFT) calculations at the PBE0-D3/def2-TZVP level of theory to investigate the NCIs present in the catalytic site. The results were analyzed using the Quantum Theory of Atoms in Molecules (QTAIM) and the Non-Covalent Interaction (NCIplot) visual index to identify and quantify S lp–π, CH–π and metal-associated (e.g., Zn SpBs) interactions present in those systems.

We believe the results reported herein (i) will provide new insights into the NCIs responsible for the stabilization of the key MET residue in Zn-SVMPs, which might have great impact in the fields of chemical biology and rational drug design and (ii) will assist in increasing the visibility of Zn-SpBs among the bioinorganic chemistry community.

Computational methods

PDB survey criteria

The PDB was inspected on 07/09/2025 to investigate the NCIs involving the Zn(II) center in SVMPs. To achieve this, the following criteria were used:

• Only X-ray structures were considered.

• No specific X-ray structure resolution criteria were imposed, although most of the structures found were solved at ∼1.0–2.2 Å resolution.

• The total number of structures found was 25.

All the retrieved structures were manually inspected and structurally analyzed to identify/analyze the presence of NCIs involving the MET-turn and the Zn(II) center, resulting in the data shown in Table 1 (for some selected examples, see below) and Table S1 (for the rest of structures, see SI). During this phase, we inspected the interactions involving the S–CH₃ group and the Zn(II) center using the following distance criteria, in line with Álvarez's van der Waals (vdW) radii definition:³⁵

• To characterize the interaction as a S lp–π bond involving the Zn coordinated HIS residues, the distance criteria used was 1.78 Å ≤ d ≤ 3.66 Å + 0.5 Å, which considers both S and C covalent and vdW radii and 1.76 Å ≤ d ≤ 3.55 Å + 0.5 Å, which takes into account both S and N covalent/vdW radii.

• To characterize the interaction as a CH–π bond involving the S–CH₃ and HIS groups, the distance criteria used was 1.46 Å ≤ d ≤ 3.54 Å + 1 Å, which considers the covalent and vdW radii of the C atom and 1.44 Å ≤ d ≤ 3.43 Å + 1 Å, which takes into account both C and N covalent/vdW radii.

• To characterize the interaction as a CH⋯Zn bond, the distance criteria used was 1.98 Å ≤ d ≤ 4.16 Å + 1 Å, which into account both C and Zn covalent/vdW radii.

• The 1 Å added in both CH–π and CH⋯Zn distance range corresponds to the typical length of a covalent C–H bond.

Table 1 List of X-ray crystal structures found during the PDB survey including the PDBID, the resolution (in Å), the type of metalloproteinase, the organism as well as S lp–π, CH–π and CH⋯Zn distances (d_S⋯HIS, d_CMET⋯HIS and d_C⋯Zn, in Å)

PDBID	Res.	Organism	dS⋯HIS ^a	dCMET⋯HIS ^a	dCMET⋯Zn ^a
a Values given as the shortest distance between the two counterparts.
PI SVMPs
1HTD	2.10	Crotalus atrox	3.69	3.67/3.47	4.44
1BUD	1.90	Deinagkistrodon acutus	3.68	3.68/3.38	4.38
1BSW	1.95	Deinagkistrodon acutus	3.64	3.97/3.36	4.78
1ATL	1.80	Crotalus atrox	3.71	3.76/3.30	4.64
1IAG	2.00	Crotalus adamanteus	3.63	3.72/3.49	4.69
4Q1L	1.90	Bothrops leucurus	3.67	3.72/3.51	4.52
3GBO	1.77	Bothrops moojeni	3.84	3.81/3.60	4.67
4J4M	1.80	Protobothrops mucrosquamatus	3.92	3.98/3.67	4.74
2W14	1.08	Bothrops asper	3.84	3.98/3.74	4.79
1YP1	1.90	Deinagkistrodon acutus	3.85	3.28/3.78	4.84
1ND1	1.93	Bothrops asper	3.90	4.00/3.69	5.02
1QUA	2.20	Deinagkistrodon acutus	3.94	3.74/3.48	4.65
1WNI	2.20	Trimeresurus flavoviridis	4.09	3.65/3.33	4.90
2W12	1.46	Bothrops asper	3.83	3.89/3.74	4.73
2W13	1.14	Bothrops asper	3.86	3.92/3.74	4.79
2W15	1.05	Bothrops asper	3.83	3.89/3.74	4.79
PIII SVMPs
3DSL	2.70	Bothrops jararaca	3.66	3.40/3.47	4.31
2DW2	2.70	Crotalus atrox	3.48	3.58/3.86	4.72
2DW1	2.50	Crotalus atrox	3.72	3.83/3.60	4.80
3HDB	2.31	Deinagkistrodon acutus	3.78	4.30/3.56	5.15
2DW0	2.15	Crotalus atrox	3.72	3.94/3.69	4.74
2ERP	2.95	Crotalus atrox	3.40	3.65/3.23	5.03
2ERO	2.50	Crotalus atrox	3.61	3.84/3.62	5.21
2ERQ	2.50	Crotalus atrox	3.96	4.13/3.45	5.14
3K7N	2.30	Naja atra	3.83	3.76/3.64	4.90

---

Modeling PDB structures and cluster model calculations

From the 25 structures retrieved from the PDB, eight of them were chosen for further phases of this study. Four of them belong to the PI family (1BUD, 1IAG, 4Q1L and 1WNI) while the rest belong to the PIII family (2DW2, 2ERO, 3HDB and 3K7N). These were selected as representative examples based on their known hemorrhagic/non-hemorrhagic activity, their respective SVMP family as well as their X-ray structure resolution.

For the selected PI and PIII SVMPs, the following computational modeling procedure was followed:

• Firstly, a spherical cutoff of 5 Å around the Zn(II) center (Zn and HIS coordinated ligands) was used to select any surrounding protein residues out from Zn the coordination sphere by means of the Discovery Studio Visualizer program,³⁶ obtaining a cluster model for each structure. This led to the obtention of the geometries containing the Zn(II) coordination complex and the main residues of the MET-turn (ILE-MET-SER).

• Secondly, the cluster models were refined using the Gaussview 5.0 software³⁷ to (i) adjust the protonation state of the system, (ii) cap the N-terminal/C-terminal regions of the protein backbone with methyl groups.

• Thirdly, these cluster model geometries were used as an input for the Semiempirical Tight Binding (xTB) calculation package,³⁸ and the protein residues relaxed while keeping the protein Cαs frozen. This approach was used to preserve the protein backbone conformation and the active site architecture. During these calculations, implicit water solvation was used by means of the Generalized Born Surface Area (GBSA) model,³⁹ implemented within the xTB code (see SI for their cartesian coordinates).

As a result of these calculations, we were able to estimate the interaction energies between the Zn(II) coordination complex and the main MET-turn residues using the supermolecule approach (ΔE_MET-turn = E_{Zn(II)center+MET-turn} − E_Zn(II)center − E_MET-turn), see Fig. 2 and Table 2.

Fig. 2 Graphical representation of the calculation process regarding ΔE_MET-turn and ΔE_MET values.

Table 2 Interaction energies between (i) MET-turn and (ii) MET residue and the Zn(II) center (ΔE_MET-turn and ΔE_MET, in kcal mol⁻¹) using xTB method and DFT calculations at the PBE0-D3/def2-TZVP level of theory, respectively. The distances between the S atom and HIS residues (d_S⋯HIS), between the –CH₃ group and the HIS residues (d_CMET⋯HIS) and between the –CH₃ group and the Zn atom (d_CMET⋯Zn) are also included in Å

PDBID	ΔEMET-turn	ΔEMET	dS⋯HIS ^a	dCMET⋯HIS ^a	dCMET⋯Zn ^a
a Values given as the shortest distance between the two counterparts.
1BUD	−22.3	−9.0	3.598	3.710/3.473	4.556
1IAG	−30.3	−9.0	3.318	3.476/3.228	4.468
4Q1L	−34.0	−8.4	3.289	3.867/3.295	4.512
1WNI	−23.7	−10.4	3.512	3.568/3.371	5.497
2DW2	−29.7	−10.7	3.438	3.501/3.239	4.190
3HDB	−37.8	−7.9	3.290	3.511/3.633	4.748
2ERO	−30.0	−9.0	3.503	3.533/3.335	5.342
3K7N	−28.5	−9.6	3.296	3.486/3.281	5.028

---

DFT calculations on reduced models

Once the relaxed geometries were obtained by means of the xTB program, reduced models containing the MET residue and the Zn(II) center were built to evaluate MET⋯Zn interaction energy (ΔE_MET values in Table 2, see SI for their cartesian coordinates).

The interaction energies were computed at the PBE0⁴⁰-D3⁴¹/def2-TZVP⁴² level of theory using the program TURBOMOLE version 7.7.⁴³ The Conductor like Screening Model (COSMO-RS)⁴⁴ was utilized within the TURBOMOLE 7.7. framework during the calculations, using water as a solvent.

In this regard, the interaction energies were calculated also using the supermolecule approximation (ΔE_MET = E_{Zn(II)center+MET} − E_Zn(II)center − E_MET), see also Fig. 2 and Table 2.

The structural view of the selected PDB structures as well as the structural comparison between the initial and final optimization states of the cluster models were carried out using the PYMOL 3.1 software.⁴⁵

The MEP surfaces involving the Zn(II) centers were also computed using the hybrid PBE0 functional and the def2-TZVP basis set by means of the Gaussian 16 software⁴⁶ and analyzed using the Gaussview 5.0 program.³⁷

The calculations for the wavefunction analysis were carried out at the same level of theory using the Multiwfn software.⁴⁷ Lastly, the NCIplot⁴⁸ isosurfaces correspond to both favorable and unfavorable interactions, as differentiated by the sign of the second-density Hessian eigenvalue and defined by the isosurface color. The color scheme is a red–yellow–green–blue scale, with red for repulsive (ρ_cut⁺) and blue for attractive (ρ_cut⁻) NCI interaction density. Yellow and green surfaces correspond to weak repulsive and weak attractive interactions, respectively. The surfaces were visualized using the Visual Molecular Dynamics (VMD) software.⁴⁹

Results and discussion

Results from the PBD inspection

The initial PDB survey identified 25 SVMP catalytic domains that fulfilled the selection criteria described in the Methods, providing a structurally diverse data set to analyze the noncovalent chemistry involving the Zn(II) center and the MET-turn (see Table 1 below). A manual inspection of these structures revealed that the conserved MET residue consistently interacts within the Zn(II)–HIS pocket, establishing several NCIs that involve (i) S lp–π contacts with the HIS imidazole rings, (ii) CH–π interactions between the MET –CH₃ group and HIS residues, and (iii) longer-range CH⋯Zn contacts that may be interpreted as SpB-like interactions. These three interacting motifs are schematically represented in Fig. 1 (see above) and are present, with differences in geometry and directionality, across both P-I and P-III SVMPs.

To explore these interactions in more detail, eight representative structures were selected for further analysis: four P-I SVMPs (1BUD,²⁸ 1IAG,¹⁴ 4Q1L, and 1WNI⁵⁰) and four P-III SVMPs (2DW2,⁵¹ 2ERO,⁵² 3HDB,⁵³ and 3K7N³³). This subset spans both highly hemorrhagic and non-hemorrhagic enzymes and includes single-domain catalytic SVMPs as well as multi-domain ectodomains. For instance, 1BUD corresponds to acutolysin A, a canonical hemorrhagic P-I toxin from Agkistrodon/Deinagkistrodon acutus venom, whereas 1IAG (adamalysin II from Crotalus adamanteus) encompasses the isolated catalytic domain of a metalloproteinase found as part of larger mosaic hemorrhagic proteins and widely used as a structural prototype for metzincins. In contrast, 4Q1L (leucurolysin-A) and 1WNI (trimerelysin II/H2-proteinase) are non-hemorrhagic P-I SVMPs: leucurolysin-A is a non-hemorrhagic fibrinogenolytic enzyme from Bothrops leucurus venom, while trimerelysin II is a major non-hemorrhagic metalloproteinase from Protobothrops (Trimeresurus) flavoviridis venom.

On the other hand, the P-III group comprises full or near-full ectodomains bearing metalloproteinase, disintegrin-like and cysteine-rich domains, including catrocollastatin/VAP2B (2DW2) and VAP1 (2ERO) from Crotalus atrox, which are potent vasculotoxic SVMPs associated with endothelial damage and hemorrhagic microvascular disruption, and AaHIV (3HDB), a strongly hemorrhagic P-IIIa SVMP from Agkistrodon acutus, together with the K-like P-III SVMP from Naja atra (3K7N), which shows distinct substrate specificity and cell-modulatory activity among elapid SVMPs.

Partial views of the catalytic centers of 1BUD, 4Q1L, and 2DW2 (see Fig. 3 below) show that in all cases the MET side chain is oriented toward the Zn–HIS pocket, with the –S–CH₃ group positioned directly above the plane defined by the three Zn-coordinated HIS residues. More in detail, in the hemorrhagic P-I toxin acutolysin A (1BUD), the S atom is located at 3.68 Å from the nearest HIS ring, while the MET –CH₃ group established two CH–π contacts at 3.68 and 3.38 Å, and the C⋯Zn distance resulted 4.38 Å. On the other hand, in the non-hemorrhagic P-I enzyme leucurolysin-A (4Q1L), the MET residue adopted a very similar orientation, with S⋯HIS and C_MET⋯HIS distances of 3.67 Å and 3.72/3.51 Å, respectively, and a slightly longer C⋯Zn contact of 4.52 Å. Finally, in the hemorrhagic P-III ectodomain 2DW2, the S⋯HIS distance was slightly shorter (3.48 Å), the C_MET⋯HIS distances observed were 3.58/3.86 Å, and the C⋯Zn distance increased to 4.72 Å. Altogether, these structures illustrate a conserved noncovalent motif in which the MET S atom and methyl group are positioned at typical S lp–π and CH–π contact distances with the HIS imidazole rings (∼3.4–3.9 Å), while the MET –CH₃ group maintains a “second-shell” CH⋯Zn contact in the 4.3–4.8 Å range.^54–56

Fig. 3 Partial views of the X-ray structures corresponding to PDBIDs (a) 1BUD, (b) 4Q1L and (c) 2DW2. The NCIs involving the MET residue and the Zn(II) center are magnified inside the square parts of the figure with indication of the S⋯HIS, C_MET⋯HIS and C⋯Zn distance values.

If this comparison is extended to the full set of SVMP catalytic domains summarized in Table 1, the same picture emerges. Across all 25 structures, S⋯HIS distances for both hemorrhagic and non-hemorrhagic enzymes are comprised within a relatively narrow interval (roughly 3.3–3.9 Å), and the associated C_MET⋯HIS and C⋯Zn distances are quite similar between the two functional groups. Some of the shortest S⋯HIS contacts in the data set are found in non-hemorrhagic enzymes, while several hemorrhagic SVMPs display S⋯HIS and C_MET⋯HIS values that are very similar from those of non-hemorrhagic counterparts. Given the geometrical similarities between the hemorrhagic/non-hemorrhagic SVMPs, it can be deduced that the MET residue is part of a robust NCI recognition pattern over the Zn(II) center whose fine tuning contributes to the modulation of catalytic behavior and hemorrhagic potential.

Geometric and energetic study on selected PI/PIII SVMPs

Fig. 4 and Table 2 gather the results from the geometry optimizations on 1BUD, 4Q1L and 2DW2 systems (see Fig. S3 and S4 in SI for the rest of the selected structures). As a result, the S⋯HIS distances lied within 3.289 to 3.598 Å, while the C⋯HIS values involving the MET –CH₃ group spanned from 3.228 to 3.867 Å (see Table 2 below). The observed distances were longer than the sums of the corresponding covalent radii but fall within, or slightly beyond, the sums of vdW radii, in line with the presence of weak yet directional S⋯π and CH–π contacts. Likewise, the C⋯Zn distances lied between 4.190 and 5.497 Å, in agreement with “second-shell” CH⋯Zn contacts rather than classical coordination bonds, thus resembling Zn SpB-like geometries.^27,28,57 Interestingly, this NCI pattern is conserved across both P-I and P-III SVMPs (see also Fig. S1 and S2 in SI for the rest of selected examples), indicating that the MET-turn is not merely a passive structural element but actively participates in the recognition pattern of the Zn–HIS pocket.

Fig. 4 Structural superposition between the initial (colored in blue) and final (colored in orange) xTB geometries in 1BUD, 2DW2 and 4Q1L structures. The protein residues/groups that changed their position more noticeably during calculations were highlighted using black arrows.

The fact that this motif persists in X-ray structures with different resolution, origin, and hemorrhagic activity strengthens the idea that subtle rearrangements in the geometry and balance of these S lp–π, CH–π, and CH⋯Zn contacts may be critical for tuning catalytic efficiency and hemorrhagic potential, in agreement with prior structure–function analyses of SVMPs and other metzincins.^7,9,10,58

In Fig. 4, the structural superpositions between the X-ray crystal structures (in blue) and the xTB-optimized geometries (in orange) are shown for 1BUD, 4Q1L and 2DW2 structures (see Fig. S3 in SI for the rest of systems). After xTB relaxation, the overall architecture of the Zn–HIS pocket and the MET-turn remained very similar to the crystallographic arrangement, with backbone deviations confined to small displacements of the loop and terminal capping groups. In the eight considered systems, the Zn(II)–HIS coordination sphere is essentially preserved, and the MET-turn backbone traces almost perfectly onto the experimental structures (as expected, since the protein Cαs were frozen as indicated in the Methods section). The residues highlighted with arrows in Fig. 4 correspond mainly to side chains located in the immediate vicinity of the MET residue (such as ILE/SER and neighboring aromatic residues) and to segments of the MET-turn loop that pivot slightly toward or away from the Zn(II) center to optimize their packing. These local rearrangements, although seem subtle, are sufficient to slightly shorten or reorient key S lp–π and CH–π contacts between MET and the Zn-coordinated HIS residues, and in some cases to adjust the geometry of the CH⋯Zn second-shell contact.

An analysis of the rest of selected examples (Fig. S3 in SI) revealed a similar picture, where the capping groups corresponding to ASP/GLU/HIS residues vary their position more noticeably. In addition, regarding the coordination environment of the Zn(II) center was also preserved in all structures, with slight changes in the water position (1IAG, 1WNI and 2ERO structures) and the catalytic GLU residue (1IAG, 3HDB). Lastly, in 3HDB a more noticeable rearrangement of both the coordinating and non-coordinating residues was observed while in 3K7N changes were observed in the dihedral angles involving the three coordinated HIS residues.

To move beyond purely structural descriptors, we evaluated the interaction energies of the NCIs established between the MET-turn main residues (ILE-MET-SER) and the Zn(II) center using both semiempirical and DFT methodologies (see Table 2 below).

As it can be observed, all the complexes exhibited negative ΔE_MET-turn and ΔE_MET values, thus indicating that the MET-turn favorably contributes to the stabilization of the catalytic Zn(II) site. Secondly, the MET-turn⋯Zn(II) interaction energies spanned from −37.8 kcal mol⁻¹ (3HDB) to −22.3 kcal mol⁻¹ (1BUD), with intermediate values around −30 kcal mol⁻¹ for most of the enzymes studied (1IAG, 2DW2, 2ERO, 3K7N) and somewhat more favorable values for 4Q1L (−34.0 kcal mol⁻¹) and 3HDB (−37.8 kcal mol⁻¹) structures. Lastly, the MET⋯Zn(II) complexes were also favorable, with moderately strong interaction energy values ranging within a narrower window (−7.9 to −10.7 kcal mol⁻¹), with most values comprised between −9.0 to −10.7 kcal mol⁻¹.

Although it should be taken as a rough estimation, from these results it can be deduced that the MET residue accounts for a sizeable fraction of the total MET-turn⋯Zn(II) stabilization, while the remaining contribution arises from the rest of the loop residues considered (ILE and SER). In fact, a simple comparison of ΔE_MET and ΔE_MET-turn values shows that the energetic contribution of MET spans from 20 to 40% across the series. For instance, in 1BUD and 1WNI, the MET residue contributes roughly 40 to 44% of the overall MET-turn⋯Zn(II) interaction energy, whereas in 4Q1L and 3HDB this fraction decreased to about 25–21%. Thus, in some enzymes the MET side chain is the principal contributor to the anchoring of the loop, while in others the surrounding residues and backbone seemed to play an increasingly important role. These results are consistent with previous mutagenesis works on metzincins, where replacing the canonical MET perturbed but did not completely abolish catalytic activity, indicating that the loop environment also participated in the stabilization of the metal center.

The geometrical data discussed above provides a structural rationale for these energetic results. In general, SVMPs with more stabilizing MET-turn⋯Zn(II) interaction energies (for example, 4Q1L and 3HDB, with ΔE_MET-turn = −34.0 and −37.8 kcal mol⁻¹, respectively) displayed S⋯HIS contacts in the lower half of the observed range (≈3.29–3.40 Å) and, in several cases, at least one relatively short C⋯HIS distance, consistent with a strengthening of S lp–π and CH–π interactions. Conversely, systems with less stabilizing ΔE_MET-turn values, such as 1BUD and 1WNI (−22.3 and −23.7 kcal mol⁻¹, respectively), tended to combine slightly longer S⋯HIS contacts (approaching ∼3.60 Å) and longer CH–π distances.

On the contrary, the C⋯Zn distances (encompassed between 4.384 and 5.497 Å) showed a weak and non-systematic relationship with the ΔE values, indicating that the CH⋯Zn (SpB-like) interaction complements, rather than dominates, the overall MET stabilization. A likely explanation could be to consider the Zn σ-hole as an electrostatic “pull” for the S–CH₃ group, being the primary driving forces the S lp–π and CH–π bonds.

Molecular electrostatic potential analyses on Zn(II) coordination sites

To gain insight into the electrostatics underlying the Zn(II) centers, we computed the Molecular Electrostatic Potential (MEP) maps for the eight selected systems at the PBE0/def2-TZVP level of theory and evaluated the electrostatic potential at the Zn–HIS pocket (see Fig. 5 and S4 in SI). More in detail, the MEP values were taken over the two HIS where the MET residue is accommodated (isovalue 0.001 a.u.; see Methods), resulting in the HISV_s values gathered in Table 3. In addition, the MEP values at the σ-hole opposite to the Zn–O_WAT bond (O–ZnV_s) and the MEP minimum and maximum values (V_s,min and V_s,max) were included in Table 3.

Fig. 5 MEP values (in kcal mol⁻¹) over the Zn–HIS pocket involving the three different Zn-coordinated centers observed (see Table 1 for the rest of MEP values). Energy values at concrete points of the surface are given in kcal mol⁻¹ (0.001 a.u.). (a) Zn center from 1BUD structure, (b) Zn center from 3HDB and (c) Zn center from 4Q1L.

Table 3 MEP values of the Zn(II) centers at the PBE0/def2-TZVP level of theory measured (i) over the two HIS residues that interact with the MET residue (HISV_s, in kcal mol⁻¹) and (ii) over the Zn–O_WAT σ-hole (O–ZnV_s, in kcal mol⁻¹) using an isocontour value of 0.001 a.u. The MEP minimum and maximum values were also included (V_s,min and V_s,max, respectively, also in kcal mol⁻¹)

PDBID	HISVs	O–ZnVs	Vs,min	Vs,max
1BUD	+64.0/+55.8	+84.1	−13	+133
1IAG	+55.2/+35.8	+83.5	−25	+132
4Q1L	+3.3/+4.4	+20.1	−82	+84
1WNI	+60.8/+58.9	+86.6	−13	+142
2DW2	+64.0/+50.8	+80.3	−6	+118
3HDB	+18.8/+11.9	+8.9	−63	+72
2ERO	+57.1/+45.2	+87.8	−13	+139
3K7N	+62.7/+39.5	+80.9	−13	+140

---

Among the eight systems considered, three different Zn(II) coordination environments were observed: (i) a distorted tetrahedral center composed by three HIS residues and one water molecule (1BUD (Fig. 5a), 1IAG, 1WNI, 2DW2, 2ERO, and 3K7N structures in SI), (ii) a Zn(II) center formed by three coordinated HIS residues and two carboxylate groups (GLU and LEU) in a trigonal bipyramid geometry (3HDB, Fig. 5b), and (iii) a Zn(II) site bound to three HIS residues and one carboxylate group in a square pyramid fashion (4Q1L, see Fig. 5c).

In the 3HIS + H₂O centers (type i), the MEP over the interacting HIS rings is predominantly positive, with HISV_s values ranging between +35 and +64 kcal mol⁻¹. Additionally, the corresponding O–ZnV_s values are also strongly positive, spanning from +80.3 to +87.8 kcal mol⁻¹, pointing to a well-defined σ-hole on Zn(II) and a strongly electrophilic Zn–HIS pocket, given the cationic nature of these Zn(II) centers.

In contrast, SVMPs in which one or two carboxylate groups directly participate in Zn(II) coordination (types ii and iii) showed a markedly different electrostatic profile. The anionic environment of 3HDB (two Zn-coordinated carboxylate ligands) attenuated the MEP value at the O–Zn σ-hole, with HISV_s values of +18.8/+11.9 kcal mol⁻¹ and an O–ZnV_s MEP value of +8.9 kcal mol⁻¹. Lastly, in 4Q1L (one carboxylate ligand), the HISV_s values dropped to +3.3/+4.4 kcal mol⁻¹ while the O–ZnV_s value was reduced to +20.1 kcal mol⁻¹. In both cases, the nearby negatively charged groups partially screen the Zn-centered σ-hole, leading to a smaller positive region in the direction of the MET residue.

Comparing the MEP data with the ΔE_MET values shown in Table 2 (see above) reveals a clear qualitative relationship between the strength of the Zn-centered MEP values and the Met⋯Zn interaction energies. More in detail, SVMPs with a strongly positive electrostatic potential at the Zn(II) site (e.g., 1BUD, 1IAG and 2DW2) tend to exhibit more favorable ΔE_MET values. In fact, 1WNI and 2DW2 combine some of the more positive V_s values with the most negative ΔE_MET values (−10.4 and −10.7 kcal mol⁻¹, respectively). In contrast, the systems in which the Zn-centered σ-hole is substantially attenuated by nearby carboxylate ligands (e.g. 4Q1L and 3HDB) show less favorable MET⋯Zn(II) interaction energies (ΔE_MET = –8.4 and −7.9 kcal mol⁻¹, respectively). Thus, while there is no strict linear correlation between V_s and ΔE_MET values, the data shown support the idea that a more positive Zn-centered potential generally enhances the electrostatic component of the MET⋯Zn(II) interaction, whereas in enzymes with a damped σ-hole, the Zn stabilizing role is increasingly complemented or compensated by S/π, CH–π interactions.

QTAIM and NCIplot studies

To further characterize the NCIs present in the systems studied herein, we analyzed the reduced MET⋯Zn(II) complexes using QTAIM and NCIplot methodologies. More in detail, the results involving three representative structures (1BUD, 4Q1L, and 2DW2) are shown in Table 4 and Fig. 6, while the corresponding analyses for the remaining SVMPs are collected in Fig. S5 and Table S1 (see SI).

Fig. 6 NCIplot analysis and QTAIM distribution of intermolecular bcps (in red spheres) and bond paths in three different MET⋯Zn complexes: (a) PDBID-1BUD, (b) PDBID-4Q1L and (c) PDBID-2DW2. The name of the NCIs involved in the MET⋯Zn recognition is also indicated, including the values of the density at the bcps that characterize those interactions in a.u. (into parenthesis). NCIplot surfaces only include intermolecular contacts between MET residue and the Zn(II) centre. NCIplot colour range −0.04 au ≤ (sign λ₂)ρ ≤ +0.04 au. Isosurface value RDG = 0.5 and ρ cutoff 0.05 a.u.

Table 4 Values of the electron density (ρ), its laplacian (∇²ρ), potential (V) and kinetic (G) energies in a.u. and the –G/V ratio at the PBE0-D3/def2-TZVP level of theory characterizing the S lp–π and CH–π interactions present in structures 1BUD, 4Q1L and 2DW2

PDBID	Interaction	ρ	∇^2ρ	V	G	–G/V
1BUD	Lp–π	0.56	1.76	−0.25	0.34	1.36
	CH–π	0.64	2.09	−0.32	0.43	1.34
	CH–π	0.48	1.49	−0.24	0.31	1.29
4Q1L	Lp–π	0.98	3.13	−0.53	0.66	1.25
	CH–π	0.71	2.46	−0.41	0.51	1.24
	π–π	1.11	3.18	−0.63	0.72	1.14
2DW2	Lp–π	0.77	2.52	−0.39	0.51	1.31
	CH–π	0.54	1.76	−0.27	0.35	1.29
	CH–π	0.66	2.22	−0.37	0.46	1.24

---

As can be observed, several intermolecular bond critical points (bcps) connect the MET residue to the Zn–HIS pocket. In all three representative systems (1BUD, 4Q1L and 2DW2), the QTAIM analysis revealed several bcps per complex, which can be assigned to specific contacts between MET and the surrounding aromatic residues. More specifically, in 1BUD (Fig. 6a) a bcp and a bond path connect the S atom of the MET side chain to the π-system of a Zn-coordinated HIS residue, corresponding to a S lp–π interaction. Additionally, two bcps and bond paths connect the C atom of the MET terminal –CH₃ group to the π-systems of two vicinal HIS residues, thus characterizing two CH–π interactions. Lastly, a bcp and a bond path connected the H atom from the MET –CH₃ group to the π-system of another neighboring HIS ring, denoting the presence of an additional CH–π bond.

In 4Q1L (Fig. 6b), the topology is slightly richer, with one bcp and a bond path connecting the H atom from the MET methyl group with a C atom from a vicinal HIS ring, denoting a CH–π contact. In addition, a bcp and a bond path connected the S atom lone pair with a N atom of a HIS ring, denoting a S lp–π interaction. Furthermore, a bcp and a bond path appeared between the aromatic rings of the Zn-bound HIS and a neighboring TRP residue, thus characterizing a π–π stacking interaction. A similar pattern was found in 2DW2 structure (Fig. 6c), with bcps and bond paths characterizing S⋯N lp–π interactions and C⋯C, C⋯N and H⋯C bcps and bond paths connecting the MET methyl group with two metal-coordinated HIS moieties, denoting the presence of three concomitant CH–π interactions.

Furthermore, in 1IAG, 1WNI, 2ERO and 3K7N a similar topological landscape was found, with several bcps that connected the MET methyl group with the Zn-coordinated HIS residues as well as a bcp and a bond path connecting the MET S atom to either a C or a N atom of a vicinal HIS residue, denoting the presence of both S lp–π interactions and CH–π bonds (see Table S1 in SI for their topological parameters). Lastly, in 3HDB structure (Fig. S5) we found the presence of three concomitant bcps and bond paths that connected the two H atoms from the MET methyl group with the three neighboring imidazole rings. Moreover, a bcp and a bond path connected the S atom with a C atom of a metal coordinated HIS.

Moreover, the Laplacian values observed were positive in all the cases, as it is common in closed shell calculations. Furthermore, the –G/V ratios extracted from QTAIM data (see Table 4) fall between 1.1 and 1.4 for all the cases considered, which is characteristic of interactions that are neither purely ionic nor strictly covalent but display an intermediate, partially shared character typical of many NCIs, such as hydrogen bonds and σ-hole interactions.^59–61

On the other hand, the NCIplot analysis provided a complementary, real-space visualization of these contacts (see also Fig. 6 above). For all three complexes, green to blue isosurfaces were observed between (i) the MET S atom and the imidazole rings (lp–π), (ii) around the MET –CH₃ group and the HIS rings (CH–π), and (iii) in the region bridging the MET –CH₃ and the Zn(II) center (CH⋯Zn). The relative size and intensity of these isosurfaces varied between 1BUD, 4Q1L, and 2DW2 structures, indicating that the balance of lp–π, CH–π, and metal-associated interactions is system dependent. In this regard, 4Q1L structure exhibited a somewhat more extensive aromatic contact surface between TRP and HIS residues, while 1BUD and 2DW2 displayed more compact but still well-defined interaction surfaces. Lastly, in the case of 3HDB structure, the NCIplot analysis was also in agreement with the results obtained from the QTAIM method and revealed the presence of long range CH⋯Zn interactions, as denoted by the greenish isosurface located between the MET methyl group and the metal ion (see also Fig. S5 in SI).

Conclusions

Herein, by analyzing high-resolution PDB structures and performing theoretical calculations on representative hemorrhagic and non-hemorrhagic SVMPs, we showed that the Met-turn is consistently oriented toward the Zn–HIS pocket and engages in an array of noncovalent interactions. In addition, through a computational at the PBE0-D3/def2-TZVP level of theory, these interactions were investigated, providing energetic and geometrical data regarding S lp–π and CH–π interactions between the MET residue and Zn coordinated HIS residues as well as CH⋯Zn contacts that exploit a Zn-centered σ-hole and can be interpreted as SpB-like interactions. These results were complemented with real-space analyses via QTAIM and NCIplot methodologies, revealing that these MET-centered contacts correspond to a network of weak but directional NCIs with clear topological and visual signatures.

Although hemorrhagic potential cannot be reduced to a single energetic/structural descriptor, our results indicate that enzymes with distinct hemorrhagic phenotypes exhibited characteristic patterns in the strength and distribution of MET-turn interactions with the Zn(II) center. We believe the results presented herein will be useful to those scientists devoted to bioinorganic chemistry and rational drug design (by designing novel small-molecule Zn(II) ligands) as well as to expand the current knowledge of SpBs among the chemical biology community.

Author contributions

C. J. R. and C. N. performed the PDB survey and the theoretical calculations while A.B. wrote the article and directed the study. All authors have given approval to the final version of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: Fig. S1 and S2: Structural view of 1IAG, 1WNI, 2ERO, 3HDB and 3K7N systems; Fig. S3: Comparison between the initial and final state of the cluster models involving 1IAG, 1WNI, 2ERO, 3HDB and 3K7N structures; Fig. S4: MEP surfaces of 1IAG, 1WNI, 2DW2, 2ERO and 3K7N systems; Fig. S5: QTAIM/NCIplot analyses of 3HDB structure; Table S1: Values of the electron density (ρ), its laplacian (∇²ρ), potential (V) and kinetic (G) energies in a.u. and the –G/V ratio at the PBE0-D3/def2-TZVP level of theory for the MET⋯Zn complexes present in structures 1IAG, 1WNI, 2ERO, 3HDB and 3K7N systems; Cartesian coordinates of the cluster models; Cartesian coordinates of MET⋯Zn models. See DOI: https://doi.org/10.1039/d6dt00716c.

This study was carried out using publicly available data from Protein Data Bank at https://www.rcsb.org/, with accession numbers: 1HTD, 1BUD, 1BSW, 1ATL, 1IAG, 4Q1L, 3GBO, 4J4M, 2W14, 1YP1, 1ND1, 1QUA, 1WNI, 2W12, 2W13, 2W15, 3DSL, 2DW2, 2DW1, 3HDB, 2DW0, 2ERP, 2ERO, 2ERQ, 3K7N.

The codes used to perform the calculations and analyze the data (including the software version) are: Discovery Studio Visualizer program version 2019 (https://www.3ds.com/products/biovia/discovery-studio/visualization), Turbomole software version 7.7 (https://www.turbomole.org/), Gaussian version 16 (https://gaussian.com/gaussian16/), MultiWFN version 2026.3.18 (https://sobereva.com/multiwfn/index.html), Gaussview version 5.0 (https://gaussian.com/glossary/g09/), NCIplot version 4.0 (https://www.lct.jussieu.fr/pagesperso/contrera/index-nci.html), Semiempirical Tight Binding (xTB) version 6.7.0. (https://xtb-docs.readthedocs.io/en/latest/xtbrelatedrefs.html), VMD version 1.9.4 (https://www.ks.uiuc.edu/Research/vmd/), PYMOL software version 3.1 (https://www.schrodinger.com/platform/products/pymol/).

Acknowledgements

A. B. thanks the MICIU/AEI for financial support (PID2023-149492NB-I00 FEDER funds) as well as CTI (UIB) for computational facilities.

References

1. D. J. Williams, M. A. Faiz, B. Abela-Ridder, S. Ainsworth, T. C. Bulfone, A. D. Nickerson, A. G. Habib, T. Junghanss, H. W. Fan, M. Turner, R. A. Harrison and D. A. Warrell, PLoS Neglected Trop. Dis., 2019, 13, e0007059.
2. S. Allen and J. K. Chakma, Trop. Med. Health, 2025, 53, 74.
3. World Health Organization, https://www.who.int/publications/i/item/9789241515641, (accessed 19 February 2026).
4. J. M. Gutiérrez, J. J. Calvete, A. G. Habib, R. A. Harrison, D. J. Williams and D. A. Warrell, Nat. Rev. Dis. Primers, 2017, 3, 17063.
5. R. Minghui, M. N. Malecela, E. Cooke and B. Abela-Ridder, Lancet Global Health, 2019, 7, E837.
6. J. W. Fox and S. M. T. Serrano, FEBS J., 2008, 275, 3016.
7. F. S. Markland and S. Swenson, Curr. Pharm. Des., 2013, 19, 4630.
8. E. Oyama and H. Takahashi, Molecules, 2017, 22, 1305.
9. C. Herrera, T. Escalante, M.-B. Voisin, A. Rucavado, D. Morazán, J. K. A. Macedo, J. J. Calvete, L. Sanz, S. Nourshargh, J. M. Gutiérrez and J. W. Fox, PLoS Neglected Trop. Dis., 2015, 9, e0003731.
10. J. M. Gutiérrez, T. Escalante, A. Rucavado and C. Herrera, Toxins, 2016, 8, 93.
11. G. F. S. Andrei, G. F. Adriana, K. Hajnal and D.-L. Muntean, Acta Med. Marisiensis, 2016, 62, 106.
12. W. Bode, F.-X. Gomis-Rüth and W. Stöckler, FEBS Lett., 1993, 331, 134.
13. F. X. Gomis-Rüth, J. Biol. Chem., 2009, 284, P15353–P13357.
14. F. X. Gomis-Rüth, L. F. Kress and W. Bode, EMBO J., 1993, 12, 4151.
15. C. Tallant, R. García-Castellanos, U. Baumann and F. X. Gomis-Rüth, J. Biol. Chem., 2010, 285, 13951.
16. T. Hege and U. Baumann, J. Mol. Biol., 2001, 314, 181.
17. A. E. Oberholzer, M. Bumann, T. Hege, S. Russo and U. Baumann, Biol. Chem., 2009, 390, 875.
18. G. S. Butler, E. M. Tam and C. M. Overall, J. Biol. Chem., 2004, 279, 15615.
19. R. Richardson-Matthews, K. Velko, B. Bhunia, S. Ghosh, J. Oktawiec, J. S. Brunzelle, V. T. Dang and A. I. Nguyen, J. Am. Chem. Soc., 2025, 147, 17433.
20. C. C. Valley, A. Cembran, J. D. Perlmutter, A. K. Lewis, N. P. Labello, J. Gao and J. N. Sachs, J. Biol. Chem., 2012, 287, P34979.
21. A. L. Ringer, A. Senenko and C. D. Sherrill, Protein Sci., 2007, 16, 2216.
22. C. A. Gibbs, D. S. Weber and J. J. Warren, Biomolecules, 2021, 12, 6.
23. M. Brandl, M. S. Weiss, A. Jabs, J. Sühnel and R. Hilgenfeld, J. Mol. Biol., 2001, 307, 357.
24. Q. Ren, W. Hao, L. Wang, K. Zhou, L. Liu and H.-C. Wu, J. Am. Chem. Soc., 2025, 147, 34118.
25. A. M. Keys, D. W. Kastner, L. L. Kiessling and H. J. Kulik, J. Biol. Chem., 2025, 301, 110379.
26. A. Bauzá, I. Alkorta, J. Elguero, T. J. Mooibroek and A. Frontera, Angew. Chem., Int. Ed., 2020, 59, 17482.
27. H. S. Biswal, A. K. Sahu, A. Frontera and A. Bauzá, J. Chem. Inf. Model., 2021, 61, 3945.
28. W. Gong, X. Zhu, S. Liu, M. Teng and L. Niu, J. Mol. Biol., 1998, 283, 657.
29. P. A. Fernandes, Toxicon, 2024, 248, 107865.
30. S. Chinnasamy, S. Nagamani and K. Muthusamy, RSC Adv., 2015, 5, 70566.
31. A. Bennacer, H. Boukhalfa-Abib and F. Laraba-Djebari, Protein J., 2021, 40, 589.
32. H. G. Wallnoefer, T. Lingott, J. M. Gutiérrez, I. Merfort and K. R. Liedl, J. Am. Chem. Soc., 2010, 132, 10330.
33. H.-H. Guan, K.-S. Goh, F. Davamani, P.-L. Wu, Y.-W. Huang, J. Jeyakanthan, W.-g. Wu and C.-J. Chen, J. Struct. Biol., 2010, 169, 294.
34. H. M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, I. N. Shindyalov and P. E. Bourne, Nucleic Acids Res., 2000, 28, 235.
35. S. Álvarez, Dalton Trans., 2013, 42, 8617–8636.
36. Dassault Systèmes BIOVIA, Discovery Studio Visualizer (version 19), San Diego, USA, 2019.
37. R. Dennington, T. A. Keith and J. M. Millam, GaussView (version 5), Semichem Inc., Shawnee Mission, KS, USA, 2016.
38. C. Bannwarth, E. Caldeweyher, S. Ehlert, A. Hansen, P. Pracht, J. Seibert, S. Spicher and S. Grimme, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2021, 11, e1493.
39. D. Qiu, P. S. Shenkin, F. P. Hollinger and W. C. Still, J. Phys. Chem. A, 1997, 101, 3005.
40. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158.
41. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
42. A. Schäfer, H. Horn and R. Ahlrichs, J. Chem. Phys., 1992, 97, 2571.
43. S. G. Balasubramani, G. P. Chen, S. Coriani, M. Diedenhofen, M. S. Frank, Y. J. Franzke, F. Furche, R. Grotjahn, M. E. Harding, C. Hättig, A. Hellweg, B. Helmich-Paris, C. Holzer, U. Huniar, M. Kaupp, A. M. Khah, S. K. Khani, T. Müller, F. Mack, B. D. Nguyen, S. M. Parker, E. Perlt, D. Rappoport, K. Reiter, S. Roy, M. Rückert, G. Schmitz, M. Sierka, E. Tapavicza, D. P. Tew, C. van Wüllen, V. K. Voora, F. Weigend, A. Wodyński and J. M. Yu, J. Chem. Phys., 2020, 152, 184107.
44. A. Schäfer, A. Klamt, D. Sattel, J. C. W. Lohrenz and F. Eckert, Phys. Chem. Chem. Phys., 2000, 2, 2187.
45. Schrödinger, LLC, The PyMOL Molecular Graphics System (version 3.1), 2025.
46. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16 (Revision A.03), Gaussian, Inc., Wallingford, CT, USA, 2016.
47. T. Lu and F. Chen, J. Comput. Chem., 2011, 33, 580.
48. 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.
49. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33.
50. T. Kumasaka, M. Yamamoto, H. Moriyama, N. Tanaka, M. Sato, Y. Katsube, Y. Yamakawa, T. Omori-Satoh, S. Iwanaga and T. Ueki, J. Biochem., 1996, 119, 49.
51. T. Igarashi, S. Araki, H. Mori and S. Takeda, FEBS Lett., 2007, 581, 2416.
52. S. Takeda, T. Igarashi, H. Mori and S. Araki, EMBO J., 2006, 25, 2388.
53. Z. Zhu, Y. Gao, Z. Zhu, Y. Yu, X. Zhang, J. Zang, M. Teng and L. Niu, Biochem. Biophys. Res. Commun., 2009, 386, 159–164.
54. M. Nishio, Y. Umezawa, J. Fantini, M. S. Weiss and P. Chakrabarti, Phys. Chem. Chem. Phys., 2014, 16, 12648.
55. V. A. Adhav and K. Saikrishnan, ACS Omega, 2023, 8, 22268.
56. S. Jena, J. Dutta, K. D. Tulsiyan, A. K. Sahu, S. S. Choudhury and H. S. Biswal, Chem. Soc. Rev., 2022, 51, 4261.
57. J. A. Pereañez, A. C. Patiño, J. M. Gutiérrez and L. M. Preciado, Med. Res. Arch., 2018, 6, 1–20.
58. C. Baldo, C. Jamora, N. Yamanouye, T. M. Zorn and A. M. Moura-da-Silva, PLoS Neglected Trop. Dis., 2010, 4, e727.
59. S. Burguera and A. Bauzá, Chem. – Asian J., 2025, 20, e00736.
60. S. Burguera and A. Bauzá, JACS Au, 2025, 5, 5438.
61. C. Nicolau, S. Burguera, M. N. Piña and A. Bauzá, Phys. Chem. Chem. Phys., 2025, 27, 9463.

---