Identification of compounds from natural Peruvian sources as potential inhibitors of SARS-CoV-2 Mpro mutations by virtual screening and computational simulations

Proteins preparation

In this study, we analyzed eight critical mutations of the SARS-CoV-2 Mpro protein, namely Y54C, N142S, T190I, A191V, S139A, G11A, R298A, and R60C. These mutations were chosen for their potential impact on the core protease’s structure and function. The Y54C mutation, identified in Malaysia, and N142S, reported in various countries, were selected because they could potentially alter the N-terminal domain’s stability and the catalytic loop’s flexibility, respectively. The T190I mutation affects the interactions within the substrate binding site, while A191V influences the dimerization dynamics, a crucial process for Mpro’s functionality. The S139A, G11A, and R298A mutations, which result in the complete loss of dimerization, are essential for Mpro’s proteolytic activity. Lastly, the R60C mutation, found in Brazil and Vietnam, affects the protein’s dynamics and the inhibitor’s ability to bind to the active site. The crystal structure of SARS-CoV-2 (PDB ID: 5RE4) reported in the Protein Data Bank (https://www.rcsb.org/ pdb/) was used. Subsequently, mutated protein sequences were prepared by replacing the amino acids at positions R298A, N142S, A191V, R60C, G11A, Y54C, T190I, and S139A. These sequences were generated by homology modelling on the SWISS-MODEL server (https://swissmodel.expasy.org) using the crystal structure of SARS-CoV-2 Mpro (PDB ID: 5RE4) as a template.

Preparation of the virtual database and screening

The search for natural products was performed at the Peruvian Natural Products Database (PeruNPDB)²⁶ online web server (first version) (https://perunpdb.com.pe/, accessed on 23 January 2022) whereas the simplified molecular-input line-entry system (SMILE) of each compound of was the upload into OpenBabel within the Python Prescription Virtual Screening Tool (PyRx)²⁷ and the subjection to energy minimization; whereas PyRx performs structure-based virtual screening by applying docking simulations using the AutoDock Vina tool.²⁸ Likewise, the FASTA sequence of the Crystal Structure of SARS-CoV-2 main protease (Mpro) (PDB: 5RE4) was subjected to a BLAST²⁹ search (accessed on 16 April 2022) whereas all the mutants were selected and subjected to automated modeling in SWISS-MODEL³⁰ server (accessed on 17 April 2022). For the analysis, the search space encompassed the whole of the modeled 3D models; and the docking simulation was then run at an exhaustiveness of eight and set to only output the lowest energy pose. Multiple sequence alignments of the Mpro and mutant sequences were visualized using the msa package (version 1.22.0)³¹ in the R programming environment (version 4.0.3). The heatmap plot was generated using GraphPad Prism version 9.4.0 for Windows, GraphPad Software, San Diego, California USA (www.graphpad.com).

Molecular dynamics simulation and Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) calculation

The simulation of the motion is realized by the numerical solution of the classical Newtonian dynamic equations. We used Gromacs v. 2020³² to calculate the molecular dynamics (MD) simulation and the AMBER-99SB-ILDN force field. The topologies for the Amber force field were determined on the ACPYPE server (https://www.bio2byte.be/acpype/) for the best metabolite against mutates Mpro. Each system was included in the centre of a cube box of 10 on each side. Likewise, water molecules were added (water model TIP4P). The energy minimization was carried out with the steep-descendent integrator with 200000 calculation steps. Herein, the MD simulation in the canonical ensemble NVT was done for a time of 1ns. Finally, the production of MD continued 100 ns in the isobaric-isothermal ensemble considering the Parrinello-Rahman barostat (1 bar) and V-rescale thermostat (309.65 K). The binding free energy estimation by MM/GBSA (Molecular Mechanics/Generalized Born Surface Area) was calculated with the suit mmpbsa.py³³ from AmberTools20³⁴ and gmx MMPBSA v1.4.1.³⁵ The equations related to calculations of binding free energies are the following:

The equation that determines the electrostatic solvation energy (∆G_GB) considers (∆E_MM) which is the variation be- tween the minimized energy of the protein-ligand complexes of the study which includes the van der Waals (∆E_vdw) and electrostatic (∆E_ele) contributions, while (∆G_SA) is the difference in surface area energies for protein and ligand and (−T∆S) refers to the contribution of entropy at temperature T.

Finally, the graphical visualizations were made with Visual Molecular Dynamics (VMD),³⁶ allowing interactive visualization with an easy-to-use interface. The interpretation of the molecular interactions was recreated with Maestro (Schrodinger) 2D interactions diagram. Likewise, the Molecular dynamics simulation results were performed by the Gromacs tools, and the values were processed by Gnuplot 5.2 (http://gnuplot.info/) command-driven interactive function plotting program.

Mutant SARS-CoV-2 Mpro description

SARS-CoV-2 Mpro is a cysteine protease of 67.6 kDa, and its structure possesses a catalytic dyad (Cys145 and His41) with a substrate-binding pocket located in a cleft between domains I and II. The secondary structure of Mpro has 10 alpha helixes, 13 beta sheets, 8 beta protrusions, 7 beta hairpins, 22 beta turns, 5 gamma turns, and 9 helix-helix interactions. In this work, we used the access code PDB ID: 5RE4 which was downloaded from the Protein Data Bank. This crystal structure was determined by the X-ray diffraction method with a resolution of 1.88 Å.

Besides, we focused on the analysis of eight mutations registered in different parts of the world. First is the Y54C mutation reported in Malaysia, and the N142S mutation was reported 17 times in 5 different countries. T190I is a mutation identified in 15 countries, such as South Africa and the USA. The mutation A191V is characterized by having an occurrence rate of 0.30% and is present in more than 34 countries. Besides, the S139A, G11A, and R298A mutation results provided a better understanding of the dimerization and catalytic mechanism of the Mpro.^24,37 In Brazil and Vietnam, the R60C mutation was reported, affecting the protein dynamics and the inhibitor’s binding within its active site.¹⁵ The R298A leads to the interruption of the dimeric conformation and irreversible inhibition of the enzyme’s catalytic activity,²⁴ and the G11A mutation avoids the insertion of the N finger region (residues 1-9) and therefore wholly declines its activity.³⁸ The location of the eight mutations is shown in Figure 1.

Figure 1. 3D representation of SARS-CoV-2 Mpro in which the eight mutations are located.

Figure 2 shows the sequence alignment of Mpro mutations. The black square selects the variation of residues by mutant Mpro. The G11A, Y54C, and R60C mutations are located close to the His41 residue and in Domain I from Mpro. Two mutations (S139A and N142S) are present in Domain II and close to Cys145, and it is expected that these protein structures could show different behaviour than Mpro without mutations. On the other hand, it was also observed that mutations in T190I and A191V are in the connection of Domain II and Domain III. For the case of R298A mutation, it can be observed near Domain III.

Figure 2. Sequence alignment of SARS-CoV-2 Mpro with the different proteins mutated.

The in silico prediction of disease-causing variants was conducted using the publicly accessible Meta-SNP tool (https://snps.biofold.org/meta-snp/).³⁹ This approach is distinguished by its integration of four established methods: PANTHER, PhD-SNP, SIFT, and SNAP each with predefined default threshold parameters. The prediction criteria are as follows: PANTHER, PhD-SNP, and Meta-SNP yield values between 0 and 1 (mutations with values >0.5 are predicted as disease-causing); SIFT produces a positive value (mutations with values >0.05 are predicted as neutral); and SNAP normalizes its output between 0 and 1 (mutations with values >0.5 are predicted as disease-causing). The predictions from Meta-SNP and its integrated tools provide a valuable initial assessment of the potential impact of nsSNVs (non-synonymous Single Nucleotide Variants). However, the variability in predictions underscores the importance of complementing these bioinformatic tools with experimental studies to achieve conclusive validation of each variant’s pathogenicity. Of particular interest are the Y54C and S139A mutations, which showed strong predictions towards pathogenicity and, therefore, warrant further in-depth analysis in future research (See Table 1).

Virtual screening analysis

The Virtual Screening technique, widely used for drug discovery, seeks to identify potential compounds for a particular therapeutic target. This approach allowed us to find new possible candidates within the PeruNPDB dataset against one of the therapeutic targets from mutant Mpro of SARS-CoV-2.

The 84 substances included in the study are taken from the original PeruNPDB collection; the most recent dataset consists of 280 substances. Figure 3 shows the gradient palette, the violet color indicated strong binding (∆G<-12 kcal/mol), while the yellow color indicated weak binding (∆G>-2 kcal/mol). In this heat map, rutin is shown as the best compounds. However, for the T190I and Y54C mutations, the color intensity is lower compared to the other mutations.

Figure 3. Heat map analysis of binding constants of metabolites from Peruvian native plants screened against mutated Mpro of SARS-CoV-2.

In Table S1 of the Supplementary Material, the values of coupling energies are reported. The values for Mpro wild, A191V, G11A, N142S, R60C, R289C, S139A, T190I, and Y54C were -10.7 kcal/mol, -10.4 kcal/mol, -10.7 kcal/mol, -10.4 kcal/mol, - 10.7 kcal/mol, -10.7 kcal/mol, -10.7 kcal/mol, -9.4 kcal/mol, and -9.1 kcal/mol, respectively. Additionally, the results of Lipinski’s rule of five⁴⁰ and ADMET (Absorption, Distribution, Metabolism, Excretion and Toxicity) prediction obtained from http://www.scfbio-iitd.res.in/ ADMETlab v 3.0⁴¹ of the 84 compounds are shown in Table S2 and Table S3 from the Supplementary Material. The Lipinski’s rule of five analysis for the majority of the compounds revealed no significant violations in terms of molecular weight, hydrogen bond donors, and acceptors. However, some compounds did not comply with the logP values. Further experimental determination of logP values may provide a better understanding of these discrepancies. The data obtained from the ADMET analysis were found to be comprehensive and informative, allowing for a better understanding of the properties to be considered in the present study.

Molecular dynamics simulations and estimation of binding free energy

The results obtained from virtual screening helped us to consider rutin as a ligand against the different Mpro mutations. Molecular dynamics simulations allow us to understand the behavior of different mutated Mpro at an atomistic level. After analyzing 100 ns of production dynamics, the convergence of each protein is observed by Root-mean squared deviation (RMSD) analysis (See Figure 2A). This result shows us that the different types of mutations achieved equilibrium; likewise, an average RMSD between 0.1 and 0.2 nm is appreciated, an acceptable value in this structural model. The Root-mean squared fluctuation (RMSF) calculates the flexibility of individual residues that make up the Mpro protein during a simulation trajectory. The RMSF per residue diagram structurally indicates which amino acids in a protein contribute the most to a molecular motion. Figure 4B highlights the area of His41 and Cys145 amino acids where the most significant fluctuation in the His41 area occurs with the R298A mutation, while the most significant fluctuation in the Cys145 area occurred in the R60C, Y54C, R298A, and N142S mutation.

Figure 4. RMSD and RMSF plots.

A.) RMSD of eight SARS-CoV-2 Mpro with rutin. B.) RMSF of the last 5 ns per residue of each SARS-CoV-2 Mpro mutated, highlighting the principal residues of the catalytic dyad (His41 and Cys145). RMSD, Root-mean square deviation; RMSF, Root-mean square fluctuation; SARS-CoV-2, severe acute respiratory syndrome coronavirus 2; Mpro, main protease.

On the other hand, Table 2 shows us the quantitative values of the RMSD, where the G11A and R60C mutations showed the lowest average RMSD value. In contrast, the average value for the Y54C mutation was higher than the others. Regarding the average RMSF values of the last 5 ns, for the 5RE4, A191V, G11A, N142S, R60C, S139A, and T190I systems, the average RMSF of the Mpro structures oscillated by 0.8 nm, while for R298A the RMSF average results in 0.09 nm and Y54C it was 0.07 nm.

Besides virtual screening and molecular dynamics simulation studies, molecular mechanics/generalized Born surface area (MM/GBSA) was performed with all frames of the MD. Table 3 indicates the average free energy values for each system. The values show a high coupling energy estimate, indicating that the interaction was carried out correctly.

The mutation R60C showed the best interaction energy (-45.09 kcal/mol) against the different systems studied. The energy values for G11A and A191V were -41.17 kcal/mol and -40.71 kcal/mol, respectively. While the systems that showed low binding energy were mutations R298A and S139A, with average values of -24.11 kcal/mol and -25.84 kcal/mol, respectively.

Likewise, the most significant energy contribution was given by the Van der Waals energies (∆E_{V DWAALS}) in the wild system, A191V, G11A, N142S, R298A, S139A, and T190I. These types of energy are weak and have short-range inter- actions; in biological systems, they play a significant role in stabilizing protein-small molecules. On the other hand, in the R60C and Y54C systems, the energy contribution is given by electrostatic energies (∆E_ELE). Electrostatic energy considers into account the charges of each atom in the system, which depend on the medium in which they are found; these have greater scope, and the force of interaction it possesses is linked to the relative orientations it accepts.

Figure 5 shows the last frame of each simulation of the Mpro-rutin complex. In general, it was observed that the interactions in the active site are due to the formation of hydrogen bonds. However, we observed some changes in the region around the active site for the mutations occurring in N142S and Y54C. The residues around N142S are mostly hydrophobic (green contour), hence N142S exhibits a greater energy contribution from hydrophobic interactions (∆E_{V DWAALS} = -54.18 kcal/mol, higher than the other mutations). In Y54C, the residues around rutin are polar (sky blue contour), demonstrating its high energetic contribution by electrostatic interactions (∆E_ELE = -66.78 kcal/mol more elevated than the other mutations).

Figure 5. 2D interaction diagram of rutin.

The pink arrow lines represent the hydrogen bond. Mpro, main protease.

Ethics and consent

Ethical approval and consent were not required.

Underlying data

Figshare: Dataset of results from virtual screening and molecular dynamics simulations. https://doi.org/10.6084/m9.figshare.27156276.v1.⁶⁵

This project contains the following underlying data:

Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).