Unraveling the Reaction Mechanism of Russell’s Viper Venom Factor X Activator: A Paradigm for the Reactivity of Zinc Metalloproteinases?

Snake venom metalloproteinases (SVMPs) are important drug targets against snakebite envenoming, the neglected tropical disease with the highest mortality worldwide. Here, we focus on Russell’s viper (Daboia russelii), one of the “big four” snakes of the Indian subcontinent that, together, are responsible for ca. 50,000 fatalities annually. The “Russell’s viper venom factor X activator” (RVV-X), a highly toxic metalloproteinase, activates the blood coagulation factor X (FX), leading to the prey’s abnormal blood clotting and death. Given its tremendous public health impact, the WHO recognized an urgent need to develop efficient, heat-stable, and affordable-for-all small-molecule inhibitors, for which a deep understanding of the mechanisms of action of snake’s principal toxins is fundamental. In this study, we determine the catalytic mechanism of RVV-X by using a density functional theory/molecular mechanics (DFT:MM) methodology to calculate its free energy profile. The results showed that the catalytic process takes place via two steps. The first step involves a nucleophilic attack by an in situ generated hydroxide ion on the substrate carbonyl, yielding an activation barrier of 17.7 kcal·mol–1, while the second step corresponds to protonation of the peptide nitrogen and peptide bond cleavage with an energy barrier of 23.1 kcal·mol–1. Our study shows a unique role played by Zn2+ in catalysis by lowering the pKa of the Zn2+-bound water molecule, enough to permit the swift formation of the hydroxide nucleophile through barrierless deprotonation by the formally much less basic Glu140. Without the Zn2+ cofactor, this step would be rate-limiting.

: Schematic representation of: (A) Topological overview of RVV-X highlighting the catalytic motif (orange shadow), the cysteine residues (yellow spheres) and the binding residues (blue spheres), (B) The FX segment docked on the RVV-X. Created with BioRender.com.   Scheme S2: Schematic representation of the final mechanism obtained for RVV-X consistent with the proposed "water-promoted pathway" mechanism.

Parametrization of the metalloproteinase's metal center
The bonded model was adopted to represent the interactions between the metal ion and the neighboring residues.
Protonation states were assigned based on the pKa values predicted by the H++ web server (http://biophysics.cs.vt.edu/H++) 1, 2 . However, H++ does not consider either the metal ion or the catalytic water while adding hydrogen atoms, so the protonation state of His139 was manually fixed. Besides the water molecule, the zinc ion is coordinated by three histidines, so only one nitrogen atom should be protonated, and the side chain should be neutral. For that reason, the extra "Hε2" atom was deleted, and the residue was renamed "HID" instead of "HIE". It is noteworthy mentioning that according to the general mechanism of action, the side chain of the Glu140 must be deprotonated to be capable of polarizing the catalytic water molecule and consequently receive a proton from it. The resulting PDB file was renumbered using pdb4amber.
In the next step, an MCPB.py script generated three models using the ff14SB force field, a small model used for bond and angle parameters calculation, a standard model containing the atom type information, and a large model for charge (RESP) calculation. The small model contained the zinc metal with the coordinating histidine's side chains, in which hydrogen atoms were added to the truncated bonds, whereas the bigger model possessed the smaller model with entire residues, and NME and ACE capped the truncated bonds.
Then, the B3LYP/6-31G* level of theory 3, 4 was employed to perform geometry optimization and force constant calculation on the smaller model. Finally, Merz-Kollman atomic charges were calculated for the larger model with the Restricted Electrostatic Potential (RESP) method and assigned to the corresponding atoms. These calculations were carried out with the Gaussian 09 software 5 .
Finally, the Seminario method, which uses the Cartesian Hessian matrix, derived the metal site force field parameters for AMBER, including the connection length, connection angle, dihedral angle, and RESP charges, and then the entire forming metal bond was integrated. MCPB.py then returned the .frcmod file necessary to generate the topology and parameters file during tLeap modeling.

Protein-peptide docking
The protein-peptide docking was done with the HPepDock online server. The software generated ten poses ranked accordingly to the lowest binding energy. After considering the interatomic distances and orientations needed for the reaction mechanism, two poses (poses B and C, with opposite orientations of the peptide, shown below) stood out as the most plausible. The peptide coordinated well with the Zn 2+ in pose C; the scissile Arg219-Ile220 peptide bond was correctly oriented concerning the catalytic residues. However, in structure B, the orientation of the scissile Arg219-Ile220 peptide bond was nevertheless less favorable for the reaction to occur, being more away from the catalytic glutamate. To further relax the docked structures and include conformational sampling, 5 replicas of both poses were simulated by MD during 170 ns (described in the following section). Several parameters were computed, such as the root-mean-square deviation (RMSD), the root-mean-square-fluctuation (RMSF), and critical protein-peptide interatomic distances. These MD trajectory analyses were performed using the CPPTRAJ module of the AmberTools18 package. The obtained results from MD confirmed that conformation C is far more suitable to the study of the reaction mechanism, as the substrate position was more stable and presented distances more suitable for catalysis.
Finally, we compared pose C with the pose of Ilomastat, co-crystallized with the D. siamensis RVV-X. The substrate (blue sticks) and Ilomastat (red sticks) share five hydrogen bonds with the surrounding residues of the "bulge-edge segment" (Gly106, Leu105, Asn103) and the "S1'-wall forming segment" (Pro165 and Leu 167), as shown in the following figure. Furthermore, the isoleucine side chain inserts into the hydrophobic S1'specificity pocket, almost perfectly superimposed on the Ilomastat chain (shaded in pink), while the arginine carbonyl group overlaps the Ilomastat carbonyl of the hydroxamate region (shaded in green) pointing to the Zinc ion.
Superimposition of the chosen conformation of the peptide (blue sticks) with the Ilomastat inhibitor (red sticks) and its interactions with the surrounding residues of the protein (yellow sticks). The Zn 2+ cofactor is represented as a slate blue sphere.

Energy minimization process
Firstly, we minimized all water molecules and then all the hydrogen atoms. Subsequently, all atoms from the side chains were minimized, and finally, an unconstrained minimization of all systems was conducted.

Molecular dynamics simulations
The molecular dynamics simulations started with all the atoms being subjected to a slow heating procedure (equilibration) in which the system's temperature gradually increased from 0 K to 300 K, using the Langevin thermostat. The heating was performed over 5 ns at constant volume, followed by 165 ns of production at constant pressure, 1 bar, recording the trajectories every 10 000 steps. The thermostat frequency of collisions was established at 1 ps -1 ; the integration step was 1 fs. Finally, the Particle Mesh Ewald method was used with a cut-off of 10.0 Å for the intermolecular interactions. The production run occurred at constant temperature and pressure, using the Berendsen barostat considering the isothermal-isobaric ensemble (NPT), allowing a more representative system density and, therefore, a more coherent simulation of the biological phenomena under study.
On the other hand, the same molecular dynamics process was conducted during 100 ns to retrieve starting structures for the multi-PES study. The trajectory was clustered using the four reactive distances as criteria.
Specifically, the Zn 2+ -HO -, Zn 2+ -scissile carbonyl oxygen, HO --scissile carbonyl carbon, and Gu140-amine nitrogen from the scissile bond. Cutoffs were used for these four distances. The cutoff values were 2.0 Å, 3.0 Å, 2.9 Å, and 4.0 Å respectively, and were applied after a pre-selection and geometry optimization. These conformations are the most reactive, and the ones used in the calculations were taken randomly within this subset.