Isolation and In Silico Prediction of Potential Drug-like Compounds with a New Dimeric Prenylated Quinolone Alkaloid from Zanthoxylum rhetsa (Roxb.) Root Extracts Targeted against SARS-CoV-2 (Mpro)

A new dimeric prenylated quinolone alkaloid, named 2,11-didemethoxy-vepridimerine A, was isolated from the root bark of Zanthoxylum rhetsa, together with twelve known compounds. The structure of the new compound was elucidated on the basis of spectroscopic investigations (NMR and Mass). The interaction of the isolated compounds with the main protease of SARS-CoV-2 (Mpro) was evaluated using molecular docking followed by MD simulations. The result suggests that 2,11-didemethoxy-vepridimerine A, the new compound, has the highest negative binding affinity against the Mpro with a free energy of binding of −8.5 Kcal/mol, indicating interaction with the Mpro. This interaction was further validated by 100 ns MD simulation. This implies that the isolated new compound, which can be employed as a lead compound for an Mpro-targeting drug discovery program, may be able to block the action of Mpro.


Introduction
Viral diseases initiated by the coronavirus (CoV) have become one of the major public health problems worldwide in the last two decades. COVID-19 is caused by SARS-CoV-2, which is a highly contagious novel coronavirus. The recent emergence of the deadly COVID-19 due to SARS-CoV-2 has created unprecedented pandemic situations around the globe, increasing the need for antiviral molecules to treat them [1,2].
There is an urgent clinical need to find new antiviral agents against SARS-CoV-2 as, except for the recently approved Paxlovid from Pfizer and some repurposed drugs like remdesivir and favipiravir, there is no approved direct acting antiviral drug that works against SARS-CoV-2 targets.
Natural products and herbal medicines provide a unique dimension in the drug development methods for antiviral medication. From a drug discovery point of view,

Characterization of Compounds
Compound 1, isolated as a yellowish mass, showed a blue fluorescent and quenching spot when examined under 366 nm and 254 nm UV light on a TLC plate. It produced a very light-yellow color when sprayed with vanillin in sulphuric acid reagent followed by heating for 5 min. It gave an orange red color after spraying with Dragendorff's reagent. The HRESIMS showed a peak at m/z 543.2480 (MH + ) measured in the positive ion mode and solved for the molecular formula C32H34N2O6. The 1 H NMR spectrum (400 MHz, CDCl3; Table 1, Figure S1) revealed two pairs of three adjacent aromatic protons at δ 7.05 dd (H-2), 7.17 dd (H-3), 7.63 dd (H-4) and δ 7.04 dd (H-11), 7.14 dd (H-12), 7.73 dd (H-13), two methoxy groups at δ 3.89 and 3.88 ppm, and two N-methyls at δ 3.93 and 3.88 ppm. In addition, the spectrum displayed three tertiary methyls at δ 1.55 (Me-6), 1.87 (Me-6), and 1.55 (Me-15) ppm in addition to seven aliphatic protons at δ 1.57 to 3.70 ppm (Ha-Hg) ( Table 1). The 13 C NMR spectrum (Table 1, Figure S2) showed thirty-two carbons altogether, including two carbonyl carbons at δ 163.9 and 163.0 ppm. The 1 H and 13 C NMR spectrum, jointly with the HSQC ( Figure S4) spectrum, revealed carbons assigned to six aromatic methines, three aliphatic methines, two methylenes, two methoxyls, two N-methyls, three tertiary methyls, four saturated quaternary carbons (two being oxygenated) and eight aromatic/olefinic quaternary carbons. All these data were similar to those of vepridimerine A [15,16], except for the absence of two methoxyl groups at C-2 and C-11. That is, compound 1 consists of two 8-methoxy-N-methyl quinolone units joined by a C10 moiety. The small coupling constant of 6.3 Hz between Hd and He suggested a cis relationship as in vepridimerine A (Tables 1 and 2). The COSY ( Figure S3), HSQC, and HMBC ( Figure S5) spectra allowed assignment of all the carbons and protons in the molecule. The relative stereochemistry was confirmed by a NOESY (Figures S6 and S7) experiment which showed cross-peaks between Hd & He, and also between Hb & Hc and Hb & Hg

Characterization of Compounds
Compound 1, isolated as a yellowish mass, showed a blue fluorescent and quenching spot when examined under 366 nm and 254 nm UV light on a TLC plate. It produced a very light-yellow color when sprayed with vanillin in sulphuric acid reagent followed by heating for 5 min. It gave an orange red color after spraying with Dragendorff's reagent. The HRESIMS showed a peak at m/z 543.2480 (MH + ) measured in the positive ion mode and solved for the molecular formula C 32 H 34 N 2 O 6 . The 1 H NMR spectrum (400 MHz, CDCl 3 ; Table 1, Figure S1) revealed two pairs of three adjacent aromatic protons at δ 7.05 dd (H-2), 7.17 dd (H-3), 7.63 dd (H-4) and δ 7.04 dd (H-11), 7.14 dd (H-12), 7.73 dd (H-13), two methoxy groups at δ 3.89 and 3.88 ppm, and two N-methyls at δ 3.93 and 3.88 ppm. In addition, the spectrum displayed three tertiary methyls at δ 1.55 (Me-6), 1.87 (Me-6), and 1.55 (Me-15) ppm in addition to seven aliphatic protons at δ 1.57 to 3.70 ppm (Ha-Hg) ( Table 1). The 13 C NMR spectrum (Table 1, Figure S2) showed thirty-two carbons altogether, including two carbonyl carbons at δ 163.9 and 163.0 ppm. The 1 H and 13 C NMR spectrum, jointly with the HSQC ( Figure S4) spectrum, revealed carbons assigned to six aromatic methines, three aliphatic methines, two methylenes, two methoxyls, two N-methyls, three tertiary methyls, four saturated quaternary carbons (two being oxygenated) and eight aromatic/olefinic quaternary carbons. All these data were similar to those of vepridimerine A [15,16], except for the absence of two methoxyl groups at C-2 and C-11. That is, compound 1 consists of two 8-methoxy-N-methyl quinolone units joined by a C 10 moiety. The small coupling constant of 6.3 Hz between H d and H e suggested a cis relationship as in vepridimerine A (Tables 1 and 2). The COSY ( Figure S3), HSQC, and HMBC ( Figure S5) spectra allowed assignment of all the carbons and protons in the molecule. The relative stereochemistry was confirmed by a NOESY (Figures S6 and S7 (Figure 2). On the basis of the above data, compound 1 was characterized as 2,11-didemethoxy-vepridimerine A. This is the first report of the isolation of this alkaloid from a natural source.      [19], [20].

Molecular Docking with MPro
Molecular docking was performed on thirteen compounds to study their interactions with the Mpro (Table 3; Figure 3) using AutoDock Vina. The highest negative binding affinity is observed for compound 1, indicating that it can potentially serve as a chemical scaffold to develop a potent Mpro inhibitor.       Non-covalent interactions are known to play a major role, as they are considered to drive protein-ligand interactions [24]. For example, hydrogen bond and hydrophobic interactions were found to be involved in the binding of the compounds with Mpro when the poses were estimated with AutoDock Vina (Figures 3 and 4). Among the thirteen compounds, the highest negative binding affinity was observed for compound 1 (−8.5 kcal/mol). Compound 1 formed five hydrogen bonds and six hydrophobic bonds with Mpro. In the compound 8-Mpro complex, the complex was stabilized by one hydrogen bond and seven hydrophobic bonds. The compound 4-Mpro complex formed a stable network with five hydrogen bonds, four hydrophobic bonds, and one electrostatic bond, whereas in the compound 5-Mpro complex, five hydrogen bonds, four hydrophobic bonds, and one electrostatic bond were detected (Figure 4). Similar interactions were shown by the other eight compounds with Mpro (Table 3 and Figure 5).  (Table 3 and Figure 5).

Result Obtained from Molecular Dynamics Simulation
For the promising Mpro inhibitors, we conducted molecular dynamics and explored the MD stability of each MD system. A 100 ns MD simulation for the complex of Mpro with new compound 1 was performed for 100 ns. The Apo form of the Mpro was also employed for the MD run. Compound 1 has lower RMSDs (0.4-2.4) for the alpha-carbon atoms than apo form (0.4-2.7), implying that compound 1 may be more stable in physiological conditions. In Figure 6A, RMSD values of compound 1 were slightly increased to 2.1 Å after 22 ns and showed a few smaller fluctuations at 26, 74, and 82 ns. Apart from this, the trajectories generated throughout the whole run were stable. Thus, higher structural stability in the Mpro-compound 1 complex was identified. On the other hand, the RMSD of Apo exhibited similar higher fluctuations over the whole run. The average RMSD for compound 1-Mpro was lower (1.6 Å) compared to apo-Mpro (2.1 Å) The Root mean square fluctuations (RMSF) property is widely used to capture protein dynamics conservation [25]. The highest degree of flexibility was detected for apo-Mpro, whereas a similar fluctuation (average 1.2 Å) was noticed for compound 1-Mpro. The lower fluctuation indicated a higher protein structural stability. In addition, the visual analysis of MD simulation trajectories showed that compound 1 was involved in major binding interactions with the hotspot residues (HIS163, GLU166, CYS145, GLY143,

Result Obtained from Molecular Dynamics Simulation
For the promising Mpro inhibitors, we conducted molecular dynamics and explored the MD stability of each MD system. A 100 ns MD simulation for the complex of Mpro with new compound 1 was performed for 100 ns. The Apo form of the Mpro was also employed for the MD run. Compound 1 has lower RMSDs (0.4-2.4) for the alpha-carbon atoms than apo form (0.4-2.7), implying that compound 1 may be more stable in physiological conditions. In Figure 6A, RMSD values of compound 1 were slightly increased to 2.1 Å after 22 ns and showed a few smaller fluctuations at 26, 74, and 82 ns. Apart from this, the trajectories generated throughout the whole run were stable. Thus, higher structural stability in the Mpro-compound 1 complex was identified. On the other hand, the RMSD of Apo exhibited similar higher fluctuations over the whole run. The average RMSD for compound 1-Mpro was lower (1.6 Å) compared to apo-Mpro (2.1 Å).
The Root mean square fluctuations (RMSF) property is widely used to capture protein dynamics conservation [25]. The highest degree of flexibility was detected for apo-Mpro, whereas a similar fluctuation (average 1.2 Å) was noticed for compound 1-Mpro. The lower fluctuation indicated a higher protein structural stability. In addition, the visual analysis of MD simulation trajectories showed that compound 1 was involved in major binding interactions with the hotspot residues (HIS163, GLU166, CYS145, GLY143, HIS172, PHE140, HIS41, THR25, MET49, MET165, and GLN189) of the Mpro protein ( Figure 6B).
The radius of gyration is a parameter that indicates protein structural compactness. During the simulation run, a lower degree of fluctuation with its consistency suggested the greater compactness and rigidity of the system [25]. Throughout the whole run, a similar radius of gyration was identified for compound 1-Mpro while comparing with apo-Mpro. The average radius of gyration of apo-Mpro was nearly same as the compound 1-Mpro ( Figure 6C). HIS172, PHE140, HIS41, THR25, MET49, MET165, and GLN189) of the Mpro protein (Figure 6B). The radius of gyration is a parameter that indicates protein structural compactness. During the simulation run, a lower degree of fluctuation with its consistency suggested the greater compactness and rigidity of the system [25]. Throughout the whole run, a similar radius of gyration was identified for compound 1-Mpro while comparing with apo-Mpro. The average radius of gyration of apo-Mpro was nearly same as the compound 1-Mpro ( Figure 6C). The expansion of protein volume is indicated by a higher solvent-accessible surface area (SASA) value, and a low fluctuation is anticipated all through the simulation time. The assessment of SASA indicated the highest value for apo-Mpro (average of 14,160 Å 2 ) compared tothe compound 1-Mpro complex (average of 14,000 Å 2 ) ( Figure 6D). The molecular surface area (MolSA) of the drug-protein complexes was also determined. In this exploration, Compound 1 revealed similar MolSA (average 14,027 Å 2 ) as apo-Mpro (Figure 7A). A noticeable difference was detected for dihedral angles of compound 1-Mpro (average 70,261) ( Figure 7B) compared to apo-Mpro. Hydrogen bonds play a vital role in molecular recognition and the stability of the whole protein structure. During the whole 100 ns, the formed intermolecular hydrogen bonds are composed of the trajectories. Compound 1-Mpro showed the comparable number of hydrogen bonds ( Figure 7C). The expansion of protein volume is indicated by a higher solvent-accessible surface area (SASA) value, and a low fluctuation is anticipated all through the simulation time. The assessment of SASA indicated the highest value for apo-Mpro (average of 14,160 Å 2 ) compared tothe compound 1-Mpro complex (average of 14,000 Å 2 ) ( Figure 6D). The molecular surface area (MolSA) of the drug-protein complexes was also determined. In this exploration, Compound 1 revealed similar MolSA (average 14,027 Å 2 ) as apo-Mpro ( Figure 7A). A noticeable difference was detected for dihedral angles of compound 1-Mpro (average 70,261) ( Figure 7B) compared to apo-Mpro. Hydrogen bonds play a vital role in molecular recognition and the stability of the whole protein structure. During the whole 100 ns, the formed intermolecular hydrogen bonds are composed of the trajectories. Compound 1-Mpro showed the comparable number of hydrogen bonds ( Figure 7C).

Binding Free Energy Analysis
To understand how the structural changes influence ligand binding, each snapshot was subjected to MM-PBSA calculation. The binding free energy of each snapshot containing a protein-ligand complex is shown in Figure 7D

Binding Free Energy Analysis
To understand how the structural changes influence ligand binding, each snapshot was subjected to MM-PBSA calculation. The binding free energy of each snapshot containing a protein-ligand complex is shown in Figure 7D. The average binding free energy of compound 1 was calculated as −23.74 KJ/mol. It showed greater free energy of binding at 89.9 ns and 93.5 ns, which was calculated as 69.53 KJ/mol and 74.73 KJ/mol, respectively.

Discussion
Compound 1 indicated the highest binding affinity of −8.5 kcal/mol. When docking results are considered with the other three compounds, together with hydrogen bond and hydrophobic bonds, compound 1 shows very strong binding interactions with the important amino acids of Mpro [26,27]. Compared to the others, through hydrogen bond and alkyl interactions, it also showed better and more appropriate interactions with the significant residues. It can promote more interactions than the others, contributing to a higher binding affinity. Amino acid residues, MET49, CYS145, GLU166, HIS163, and MET165 in the catalytic domain of Mpro were found to participate in non-covalent interactions with most of the compounds (Figure 4 and Table 3). By MD simulation, compound 1 showed a lower average RMSD values (1.6) than apo-Mpro, indicating more structural stability with Mpro in bound states. From the RMSF values, it was observed that compound 1 had more acceptable binding stability with Mpro ( Figure 6B), which demonstrates important interaction with significant residues and their compactness in a complex form with Mpro. The compound 1 was found to be stable ( Figure 6C,D), as revealed by the Rg and SASA values. Furthermore, the compound 1-Mpro complex showed a comparable result for MolSA ( Figure 7A). All analyses from the MD simulations support the docking results, signifying that our isolated compound 1 (2,11-didemethoxy-vepridimerine A) has potential to inhibit the activity of Mpro of SARS-CoV-2.

Discussion
Compound 1 indicated the highest binding affinity of −8.5 kcal/mol. When docking results are considered with the other three compounds, together with hydrogen bond and hydrophobic bonds, compound 1 shows very strong binding interactions with the important amino acids of Mpro [26,27]. Compared to the others, through hydrogen bond and alkyl interactions, it also showed better and more appropriate interactions with the significant residues. It can promote more interactions than the others, contributing to a higher binding affinity. Amino acid residues, MET49, CYS145, GLU166, HIS163, and MET165 in the catalytic domain of Mpro were found to participate in non-covalent interactions with most of the compounds (Figure 4 and Table 3). By MD simulation, compound 1 showed a lower average RMSD values (1.6) than apo-Mpro, indicating more structural stability with Mpro in bound states. From the RMSF values, it was observed that compound 1 had more acceptable binding stability with Mpro ( Figure 6B), which demonstrates important interaction with significant residues and their compactness in a complex form with Mpro. The compound 1 was found to be stable ( Figure 6C,D), as revealed by the Rg and SASA values. Furthermore, the compound 1-Mpro complex showed a comparable result for MolSA ( Figure 7A). All analyses from the MD simulations support the docking results, signifying that our isolated compound 1 (2,11-didemethoxy-vepridimerine A) has potential to inhibit the activity of Mpro of SARS-CoV-2.

General Experimental Procedures
In the present work, NMR spectra were measured at 400 MHz for 1 H NMR spectra and 100 MHz for 13 C on a Bruker 400 TM ASCEND spectrometers (Bruker UK, Coventry, United Kingdom). The sample was dissolved in CDCl 3 -d. Chemical shifts (δ) are reported on a ppm using tetramethyl silane (TMS) and/or residual solvent peak as an internal reference. Coupling constants (J) are given in hertz. High resolution electrospray ionization mass spectrometry (HRESIMS) was used for obtaining the mass spectrum. Column chromatography (CC) was carried out using powdered silica gel 60 H (Kieselgel 60, mesh 70-230). The glass column, having a length of 92 cm and diameter 4.5 cm, was packed with silica gel. Thin layer chromatography (TLC) was obtained using pre-coated silica gel plates (silica gel 60 F254, aluminum sheets, E-Merck, Germany) and the established chromatogram was detected under UV light (366 nm and 254 nm). The glowing and quenched spots were marked, and the recognition was then visualized by spraying with Vanillin-H 2 SO 4 and modified Dragendorff's reagent.

Plant Material
The root bark of Zanthoxylum rhetsa was collected from Dhaka, Narsingdi district of Bangladesh in August 2013. The plants were identified by an expert at the Bangladesh National Herbarium and where a voucher specimen was deposited (DACB accession number is 42,528). The plant parts were cleaned properly. The parts were cut into small pieces and subjected to shade drying for a week. The dried plant part was then crushed into coarse powder by a high-capacity grinding machine with proper care.
The entire pet ether soluble fraction (9 g) was exposed to column chromatography (CC) by means of silica gel by gradient elution with pet ether/dichloromethane/ethyl acetate/methanol [PE-CH 2 Cl 2 gradient (100:0-0:100, v/v) to CH 2 Cl 2 -EtOAc (99:1-0:100, v/v) to EtOAc-MeOH (99:1-0:100, v/v)] to afford 551 fractions each with 20 mL. The Optimization at the semi-empirical PM6 method (gaseous phase) of the structure of the selected compounds was employed [28]. Later, the vibrational frequencies were computed and the absence of imaginary frequencies confirmed the stationary points correspond to minima on the Potential Energy Surface. From the RCSB Protein Data Bank (PDB ID: 6LU7) the crystal structure of the Mpro was taken. After that it was energy-minimized through a Swiss PDB Viewer [29]. The BIOVIA discovery studio (version 4.5, Vélizy-Villacoublay, France) software package was used to eliminate all of the hetero atoms and water molecules present in the crystal structure [30]. Molecular docking between Mpro with thirteen compounds was performed by AutoDock Vina protocol for the potential interactions and activity predictions [31]. The Mpro active site was set around in the docking grid box, where the center was around X = −9.83, Y = 14.38, Z = 68.48 and the dimensions were X: 24.69, Y: 31.64, and Z: 29.97. Finally, for detecting the non-covalent interaction in the docked drug-protein complex, BIOVIA discovery studio (version 4.5, Vélizy-Villacoublay, France), were visualized and detected.

Molecular Dynamics
Along with the apo form of Mpro, MD simulations of new compound 1 was performed using the YASARA software package, 21.8.27 (Vienna Austria). These simulations were executed to study the virtual stability of the ligand locate in the binding pocket. A dynamics system consisted of one copy of Mpro, one copy of docked compound, water molecules, and 0.9% NaCl at 298K temperature. The entire system was neutralized. The MD system was first relaxed for a protein-ligand complex via a series of procedures such as the steepest descent minimization and simulated annealing minimization of the solvent. For describing the macromolecular system, the AMBER14 force field was utilized. Three-dimensional periodic boundary conditions (PBC) were applied wherever the cell size was 20 Å larger than the Mpro size in all cases. Throughout all simulations, the Particle-Mesh-Ewald (PME) method was employed for long-range electrostatic interactions (with a cut-off of 8 Å). The simulation temperature was controlled by the Berendsen thermostat process. For performing 100 ns MD simulation, a time step of 1.25 fs was maintained. The snapshots were taken at every 100 ps. Depending on previous MD data analysis process root-meansquare deviation (RMSD), root-mean-square fluctuation (RMSF), radius of gyration (Rg), and solvent accessible surface area (SASA), molecular surface area (MolSA), dihedral angle, and total number of hydrogen bonds data were retrieved from the MD simulations [32][33][34]. At that point, 100 snapshot conformers were generated and saved in PDB format at every 1 ns of 100 ns MD simulation.

Binding ENERGY Calculation MM-PBSA
By using YASARA software, Gibbs free energy was calculated for analyzing the binding affinity of ligands. Interaction energy for Mpro with compound 1 was calculated. MM-PBSA methods were implemented for 51 to 100 ns.

Conclusions
In this study we identified a new natural compound, 2,11-didemethoxy-vepridimerine A (1), the dimeric prenylated quinolone alkaloid, that showed promising interaction with SARS-CoV-2 MPro. Mpro, coronavirus' main protease, is a critical enzyme mediating the production of CoV functional proteins. The molecular docking results were supported by the MD simulation calculations. This suggests that the isolated compound 1 can potentially inhibit the activity of Mpro, and can be used as a lead compound for an Mpro targeting drug-discovery program. However, these results need to be validated using in vitro wet lab experiments to fully understand the potential utility of this chemical class.