Active site-specific quantum tunneling of hACE2 receptor to assess its complexing poses with selective bioactive compounds in co-suppressing SARS-CoV-2 influx and subsequent cardiac injury

Objective: This research aims to study the target specificity of selective bioactive compounds in complexing with the human angiotensin-converting enzyme (hACE2) receptor to impede the severe acute respiratory syndrome coronavirus 2 influx mechanism resulting in cardiac injury and depending on the receptor’s active site properties and quantum tunneling. Materials and Methods: A library of 120 phytochemical ligands was prepared, from which 5 were selected considering their absorption, distribution, metabolism, and excretion (ADMET) and quantitative structure–activity relationship (QSAR) profiles. The protein active sites and belonging quantum tunnels were defined to conduct supramolecular docking of the aforementioned ligands. The hydrogen bond formation and hydrophobic interactions between the ligand–receptor complexes were studied following the molecular docking steps. A comprehensive molecular dynamic simulation (MDS) was conducted for each of the ligand–receptor complexes to figure out the values – root mean square deviation (RMSD) (Å), root mean square fluctuation (RMSF) (Å), H-bonds, Cα, solvent accessible surface area (SASA) (Å2), molecular surface area (MolSA) (Å2), Rg (nm), and polar surface area (PSA) (Å). Finally, computational programming and algorithms were used to interpret the dynamic simulation outputs into their graphical quantitative forms. Results: ADMET and QSAR profiles revealed that the most active candidates from the library to be used were apigenin, isovitexin, piperolactam A, and quercetin as test ligands, whereas serpentine as the control. Based on the binding affinities of supramolecular docking and the parameters of molecular dynamic simulation, the strength of the test ligands can be classified as isovitexin > quercetin > piperolactam A > apigenin when complexed with the hACE2 receptor. Surprisingly, serpentine showed lower affinity (−8.6 kcal/mol) than that of isovitexin (−9.9 kcal/mol) and quercetin (−8.9 kcal/mol). The MDS analysis revealed all ligands except isovitexin having a value lower than 2.5 Ǻ. All the test ligands exhibited acceptable fluctuation ranges of RMSD (Å), RMSF (Å), H-bonds, Cα, SASA (Å2), MolSA (Å2), Rg (nm), and PSA (Å) values. Conclusion: Considering each of the parameters of molecular optimization, docking, and dynamic simulation interventions, all of the test ligands can be suggested as potential targeted drugs in blocking the hACE2 receptor.


Introduction
The global severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) pandemic has been deteriorating over time considering the frequencies of spreading and the fatality ratio. Like specific disease-causing microbes, the genetic makeup of SARS-CoV-2 also keeps changing over time, and seven spontaneous mutations have already been reported [1]. There are four types of coronaviruses which belong to the Coronaviridae family, including α-CoV, β-CoV, δ-CoV, and γ-CoV [2]. Coronavirus is a zoonotic virus which encodes many open reading frames because of a positive, single-stranded RNA genome [3]. Although SARS-COV-2 belongs to the Betacoronavirus genus, it is more dangerous than the others of the same genus [4]. According to the most recently published data, the mortality rate is higher among aged people, especially over 60 years. Coronavirus disease 2019 (COVID-19) triggers various acute respiratory distress and respiratory failure, leading to cardiac injury, heart failure, and dysfunction in multiple organs [5]. SARS-CoV-2 possesses a total of five genes in its genome, such as the ORF1ab gene that encodes 16 non-structural proteins, envelope (E) gene that encodes envelope protein, spike (S) gene that encodes spike protein, membrane gene that encodes membrane (M) protein, and nucleocapsid gene that encodes nucleocapsid protein [6]. The "spike" proteins exist in two unique conformations, namely pre-fusion and post-fusion. Activation of the S protein is required to transform into its post-fusion confirmation from the pre-fusion status, which leads to membrane fusion which guides coronavirus entry into host cells. S1 and S2 are the two subunits of the spike protein, where the S1 subunit adheres to the cell surface receptor through its receptor-binding domain (RBD). Then the virus fuses with the host membrane through the S2 subunit. Diversified host receptor recognizing capabilities of the S1 subunit is responsible for the multiple variants of SARS-CoV-2. The "S protein-human angiotensin-converting enzyme (hACE2) receptor" complex formation is the main precursor of viral entry and propagation inside the host cells [7]. Individuals with strong immunity show very mild symptoms, whereas people with hypertension, heart and autoimmune diseases, respiratory tract illness, multiple organ damage such as liver, kidney, gastrointestinal tract, and central nervous system become very susceptible to SARS-Cov-2 [8]. Following SARS-CoV-2 infection, several physiological effects become more transparent and acute, like nausea, fever, coughing, vomiting with abdominal pain, and so on [9,10].
One of the potent inhibitors of the renin-angiotensin system (RAS) is ACE2, which maintains the equilibrium of blood pressure and fluid balance [11]. The ACE2 receptors are responsible for ceasing the detrimental effects caused by angiotensin II (Ang II), such as vasoconstriction, inflammation, and fibrosis, through degradation of Ang II. After degradation of Ang II, it is converted into Ang I. Ang I acts as a vasodilator and anti-proliferator [12]. ACE2 opposes the actions of Ang II either by indirectly reducing the Ang II synthesis in tissues through cleavage of Ang I or direct hydrolyzing Ang II [13]. Also, their expression varies according to age, sex, and underlying diseases like diabetes, hypertension, and cardiovascular and pulmonary disease [14]. Trimers of the spike protein present in the host receptor ACE2 arbitrates the attachment of SARS-CoV-2 to the cell membrane, which assists the virus' entry into cells [15]. The virus' entry through the receptor in the upper respiratory system and lungs significantly damages the cardiac system. It indicates that this virus plays a pathological role in myocardial ACE2 expression [16]. Studies have shown that SARS-CoV-2 has more affinity toward ACE2 receptors, along with other reports depicting those organs are receptive to SARS-CoV-2 infection, demonstrating higher levels of ACE2 expression [17]. Attachment of the spike protein with the ACE2 receptor leads to a decrease in ACE2, which results in an upregulation of Angiotensin II, causing an imbalance in RAS, which is unfavorable for the patients. This dysregulation implicates the onset and pathogenesis of hypertension [18]. Hence, the presence of the host ACE-2 receptor is the gateway for the SARS-CoV-2 virus' entry into host cells which has been highlighted in many studies recently [1]. Thus, variations present in spike protein and host ACE-2 receptor binding sites can considerably decrease the virus' effects on the patient's physiology [19]. Ang1-7 has many diverse roles, such as anti-inflammation, antioxidant, vasodilatory, and natriuretic effects regulated by the G-protein-coupled receptor [20], which can be misregulated S-hACE2' complexing creating heart injury.
Considering all of the factors above, the current study aims to identify the interactions and target specificity than 2.5 Ǻ. All the test ligands exhibited acceptable fluctuation ranges of RMSD (Å), RMSF (Å), H-bonds, Cα, SASA (Å 2 ), MolSA (Å 2 ), Rg (nm), and PSA (Å) values. Conclusion: Considering each of the parameters of molecular optimization, docking, and dynamic simulation interventions, all of the test ligands can be suggested as potential targeted drugs in blocking the hACE2 receptor.
of selective biosynthetic ligands to the hACE2 receptor, depending on its active site detection and quantum tunnel profiles. The hydrogen bindings and hydrophobic interactions responsible for blocking the S-hACE2 receptor complex formation were also analyzed. Besides, to figure out the significance of molecular dynamic simulation in evaluating the efficacies of the test bioactive components was also conducted to impede SARS-CoV-2 influx mechanism and subsequent heart injuries, sophisticated computer programming, and biostatistical algorithms were used.

Construction of the library of phytochemical compounds
A library of phytochemical aromatic compounds was prepared with 120 ligands, followed by a thorough literature review. Ligand screening has been conducted through comprehensive profiling of their physiochemical properties based on PubChem (https://pubchem.ncbi.nlm.nih.gov/), a National Center for Biotechnological Information affiliated directory of chemical substances and biological assays.

Screening the ligand library via ADMET and QSAR profiling for ligand validation
PubChem (https://pubchem.ncbi.nlm.nih.gov/) has been utilized for collecting all the 3D structures of the targeted ligands in the form of structure data file (SDF). To study their pharmacokinetic features, such as absorption, distribution, metabolism, excretion, and toxicity, all the 120 ligands have been checked via the absorption, distribution, metabolism, and excretion (ADMET) profiling that interprets the ligand's propensity inside the body [21]. For generating the ADMET profile of the selected test ligands, "Swiss ADME" (http:// www.swissadme.ch/index.php) and "Molinspiration Cheminformatics" (https://www.molinspiration.com/cgibin/properties) was employed. For secondary identification of ADMET, "pkCSM" (http://biosig.unimelb.edu.au/pkcsm/ prediction) was implemented where the ligands' toxicity parameters were emphasized. Subsequently, the potential ligands were run via "admetSAR 2" (http://lmmd.ecust.edu. cn/admetsar2/) for determining the quantitative structureactivity relationship (QSAR). After considering the ADMET and QSAR profiling, six ligands that showed promising values and properties in terms of their ADMET and QSAR were selected for supramolecular docking and molecular dynamic simulation against the control ligand serpentine. All the ligands were optimized using University of California San Francisco (UCSF) Chimera (version 1.14) [22]. The Gasteiger approach was applied to minimize energy by lowering the accumulative charge on ligands to zero [23]. The optimized ligands were then converted into a "mol2 file" for conducting molecular docking.

Preparation of the macromolecule
The 3D crystal structure of human angiotensin-converting enzyme-related carboxypeptidase (hACE2) was gathered from the database of protein data bank (PDB ID: 1R4L and Resolution 3.00 Å). The protein's crystal structure was optimized using UCSF Chimera (version 1.14) to obtain proper orientation, size, and rotations [23][24][25]. The non-standard amino acid, ions, water molecules, and ligands were deleted from the protein crystal structure during optimization to avoid interacting with the undesired parts of the receptor. In addition, to assure congenial performance during docking, missing hydrogen atoms were added to the macromolecule [26]. The minimized energy was calculated by YASARA (https://www.yasara.org) [27].

Active site prediction of the receptor macromolecule
The supramolecular docking pose of the optimized hACE2, indicating the best active site, was determined and validated using the COACH-D (https://yanglab.nankai.edu. cn/COACH-D/) algorithm [28]. Using the algorithm, a total of three suggestive binding poses were found, among which the best one was selected, considering the number of amino acid residues involved and the projected binding affinity (kcal/mol). Finally, the best active site predicted file was subjected to quantum tunneling.

Quantum tunneling on the best active site of the receptor
To gain a better understanding of the intended modalities of ligand accommodation within the hACE2 receptor, various protein tunnels were found using the Schrodinger algorithm and the CAVER Web 1.0 tools (https://loschmidt. chemi.muni.cz/caverweb/). A total of 24 tunnels were found from the initial prediction. Finally, considering the tunnel length and radius (Å), the six most viable tunnels were determined using CAVER 3.0 [29]. The tunnels were studied to validate the pre-identified best active site of the receptor for proceeding a successful supramolecular docking because protein tunneling explains efficient electron transport via protein junctions which is a precursor of super-docking [30].

Point-specific molecular docking
Molecular docking of the selected optimized ligands was undertaken using PyRx version 0.8 to scrutinize the binding efficacy of the optimized protein-ligand complexes [31]. The desired ligands and macromolecule were transformed to the format "pdbqt during the molecular docking operation. The root mean square deviation (RMSD) (Å) and binding affinity (Kcal/mol) output files were saved as comma-separated values files.

Post-docking analysis
For the initial visualization and qualitative receptor-ligand interactions analysis, Discovery Studio Visualizer (version 3.0) and PyMOL (version 2.4.1) were used sequentially, and the output files were saved as PDB files. Afterward, the quantitative hydrophobic interactions and the number of hydrogen bond formations between each of the protein-ligand complexes were analyzed using LigPlot+ (version 2.2) [32] before conducting the molecular dynamic simulation.

Molecular dynamics simulation
In the beginning, the ligand-free macromolecule 1R4L was subjected to dynamic simulation for 10 ns to investigate its natural physical alteration and its reaction with surrounding ions and water molecules utilizing the webbased dynamic simulator CABS-flex 2.0 web-based simulation (http://biocomp.chem.uw.edu.pl/CABSflex2/) [33]. Subsequently, the protein-ligand complexes were preliminarily subjected to the ligand and receptor molecular dynamics (LARMD) simulation system (http://chemyang. ccnu.edu.cn/ccb/server/LARMD/index.php) [34]. LARMD was run for 3.1 ns to understand the outcomes of each protein-ligand complex independently in terms of solvent accessible surface area (SASA), Debye-Waller factor for thermostability (B-factor), root mean square fluctuation (RMSF), principal component analysis, and RMSD analysis [32]. Finally, molecular dynamic simulation of the ligandreceptor complexes was conducted up to 20 ns using Desmond (Desmond, Schrödinger, LLC, NY) software package to investigate RMSD (Å), RMSF (Å), Rg (nm), H-bonds, SASA (Å 2 ), molecular surface area (MolSA) (Å 2 ), and polar surface area (PSA) (Å) of each protein-ligand complex [35]. The box dimension was fixed at X:Y:Z, and nullifying ions (Na+) were added as needed to get the desired results. The probe radius was adjusted to 1.4 Å to investigate the molecular surface area and solvent-accessible surface area.

Statistical analysis and graphical representation
The protein-ligand complexes have been improved and visualized using the software packages listed above. In addition,  [36,37] and R programming (version R-4.0.2 for Linux) [38,39].

Protein active site detection and quantum tunneling
The energy level of the hACE2 receptor protein was minimized from −52407.7 to −113291.2 KJ/mol, as a result of the optimization process (Fig. 1A). The three protein active sites showed their projected binding affinities as −4.1kcal/ mol ( Fig. 2A), −4.2Kcal/mol (Fig. 2B), and −6.9Kcal/mol (Fig. 2C). Considering the number of amino acid residues and binding energy, the third pose of the protein active site ( Fig. 2C) was taken for quantum tunneling. For the third active site position, six tunnels were identified (each at 64 Å length) with the radius (Å) ranged between 0.83Å and 3.17Å (Fig. 3). The tunnels are the suggested routes for molecular super docking.

Pharmacokinetic profiles of the ligands
In this study, five ligands, including control, were picked from a library of 120 renowned pharmacophores following their QSAR and ADMET profile analysis, which are piperolactam A (Fig. 1B), apigenin (Fig. 1C), isovitexin (Fig. 1E), and quercetin ( Fig. 1F) as test ligands, whereas serpentine ( Fig. 1D) was the control (Table 1). Based on the physicochemical analysis, there were no violations of Lipinski's rules in any of the tested ligands. While observing the pharmacokinetics analysis, only piperolactam A among the ligands showed partial AMES toxicity. In addition to an excellent excretion rate, these five ligands had maximum tolerated dosages ranging from −0.219 to 0.649 log mg/kg/day. The intestinal absorption of serpentine (control), apigenin, and piperolactam A were typically > 90%, whereas quercetin and isovitexin were 77.207% and 64.729% respectively. The blood-brain barrier range of the ligands, namely apigenin, isovitexin, piperolactam A, and quercetin, were −0.734, −1.375, −0.397, and −1.098, respectively, but control serpentine was 0.257. Except for serpentine, all ligands demonstrated hepatotoxicity and their LD50 ranged between 2.471 and 3.675 (Table 1).

Supramolecular docking
Supramolecular docking showed that isovitexin exhibited the highest binding affinity (−9.9 Kcal/mol) toward the enzyme, which was also higher than the control ligand, serpentine (−8.6 Kcal/mol). Besides isovitexin, quercetin also showed a higher binding affinity value (−8.9 Kcal/ mol) than the control ligand. On the other hand, apigenin      (Table 2).

Molecular dynamic simulation (20 ns)
In total, eight parameters have been considered to check the fluctuation profiles of the docked complex from the molecular dynamic simulation (Fig. 7). During 20 ns of the simulation period, 596 interactive amino acid residues of 1R4L were harvested among 1001 different frames.
From the molecular dynamic simulation analysis, the control ligand serpentine was found to have a range of RMSD values within 0 Ǻ-2.07 Å. Isovitexin displayed the  highest RMSD value of 3.752Ǻ in comparison to all test and the control ligands. On the contrary, piperolactam A showed the lowest RMSD value at 2.057 Å. However, quercetin showed higher and apigenin exhibited lower RMSD values at 2.196 Å and 2.062 Å, respectively, than the control serpentine (Fig. 7A). Similarly, isovitexin showed a significantly higher RMSF fluctuation range from 0.520 Å to 4.268 Å, as compared to the control ligand (0.424 Å-3.704 Å) and the rest of the three experimental ligands (Fig. 7B). Besides, apigenin (0.417 Å-4.054 Å) and quercetin (0.427 Å-3.944 Å) also exhibited higher RMSF fluctuations compared to serpentine. On the other hand, no significant difference in fluctuations was observed between control and piperolactam A (0.428 Å-3.883 Å). An intramolecular hydrogen bond was absent in the docked complexes containing the control ligand, serpentine, and piperolactam A (Fig. 7C). On the contrary, quercetin and apigenin possessed only one intramolecular hydrogen bond. The isovitexin-1R4L complex exhibited the highest number of intramolecular hydrogen bond interactions (Fig. 7C).
For the control ligand serpentine, the range of interactive alpha carbon atoms (Cα) was observed between 0.402 and 3.443, which was found to be very close to the apigenin (0.396-3.517) and piperolactam A (0.412-3.579).
On the other hand, a higher range of interactive alpha carbon atoms was observed in both isovitexin-1R4L (0.502-4.573) and quercetin-1R4L (0.416-4.056) complexes, compared to the serpentine-1R4L complex. However, isovitexin showed the highest interactive alpha carbon atoms range during the simulation period with the protein (1R4L) compared to all other ligands (Fig. 7D).
In this study, SASA revealed the level of exposure of 597 amino acid residues of the enzyme 1R4L, associated with ligand to solvent (water) and pharmacokinetic strength of lead molecules, including quercetin, isovitexin, piperolactam A, and apigenin along with the control ligand, serpentine. In all cases, the water probe radius of 1.4 Å and the ratio of the total area to energy for polar region 9598.23 and apolar region 14050.71 were observed with no gradient of calculation and no unknown area (  (Fig. 7E).
To determine the MolSA, a probe radius of 1.4 Å (equivalent to the van der Waals surface area of a water molecule) When compared to all test ligands, including the control (Fig. 7F), isovitexin demonstrated the highest MolSA value (364.806 Å 2 ), whereas piperolactam A exhibited the lowest (242.419 Å 2 ). However, quercetin and apigenin displayed lower molecular surface area values of 256.595 Å 2 and 244.140 Å 2 , respectively, than that of serpentine (329.282 Å 2 ).
In the molecular dynamic simulation course, the highest fluctuation of the radius of gyration (Rg) value was observed in the docked complex containing isovitexin, varying in the range between 4.760 and 4.992 nm. On the contrary, the ligand piperolactam A exhibited the lowest score between 3.111 and 3.211 nm. The control ligand, serpentine, displayed Rg values ranging from 3.913 to 4.164 nm. Besides, quercetin (3.696-3.849 nm) and apigenin (3.597-3.743nm) generated Rg values of a lower range than that of the control ligand (Fig. 7G).
Although all ligands exhibited a higher range of polar surface area than the control, isovitexin showed the widest range (297.167 Å-356.201 Å) compared to all the test and control ligands. The PSA values found for serpentine ranging from 67.614 Å up to 86.226 Å. Notable fluctuations were also observed among quercetin, piperolactam A, and apigenin with the PSA values of 284.805 Å, 129.34 Å, and 200.676 Å, respectively (Fig. 7H).

Discussion
The World Health Organization identified the Delta variant [Phylogenetic Assignment of Named Global Outbreak lineage designation B.1.617.2], initially found in India, as a variant of concern in May 2021, and linked it to an epidemic resurgence in the United Kingdom. The lineage comprised three subtypes (B1.617.1, B.1.617.2, and B.1.617.3), each with a different set of "Spike mutations" in the N-terminal domain (NTD) and the RBD that might help them evade the immune system. "B.1.617.2," commonly known as variation Delta, is thought to spread more quickly than other variants [40].
From the understanding of the pathophysiology of the SARS-CoV-2, it can be said that the spike protein acts as a ligand by binding with the hACE2 receptor facilitating the viral particle entry into the host cell. Thus, through modulator proteins as ligands, it can be possible to allosterically regulate the spike protein's binding activity [41]. Phytochemicals are naturally synthesized in plants and are well known to have numerous therapeutic properties, having several studies reported on its competence against diseases, additionally offering a variety in compounds chosen [42]. The binding of these phytochemicals with host protein ACE2 as a non-competitive molecule can confer antiviral efficacy by disrupting spike protein binding to the hACE2 receptor [43]. For the current in silico study, the chosen pharmacophore compounds against the ACE2 receptor are quercetin, isovitexin, piperolactam A, and apigenin, keeping SERPENTINE as the control ligand ( Fig.  1B-F).

ADMET and QSAR
From the assessment of the physicochemical properties, only isovitexin among the five ligands was observed to partially violate the Lipinski rule of 5, leading to its deflection than the other ligands exhibited drug-likeness properties [25]. All ligands displayed prominent excretion rates along with acceptable ranges of maximum tolerated dosages. More than 90% intestinal absorption was witnessed for the ligands serpentine (control), apigenin, and piperolactam A, whereas lower than 90% was displayed by quercetin and isovitexin. AMES toxicity was tested positive for piperolactam A only. The blood-brain barrier range of the ligands apigenin, isovitexin, piperolactam A, and quercetin came out as negative, whereas the value for serpentine only came out as positive. With the control ligand being the exception, the hepatotoxicity came negative for the rest of the four ligands (Table 1). In the case of the LD50 dosage, serpentine acted abnormally again, whereas the LD50 range for the other four ligands were acceptable ranges.

Protein active site detection and quantum tunneling
A web server named "CASTp" (http://sts.bioe.uic.edu/ castp/index.html?3trg) was used for the prediction of the active binding site of the ACE2 receptor, but it was able to give only a qualitative-predicted result (Fig. 1A). COACH-D algorithm was utilized for the quantitative assessment of the predicted active site of the hACE2 receptor protein [28,34], which generated quantitative data of a total of eight best calculating binding poses. Upon analyzing the binding energy and amino acid residue positions, the top three binding poses were chosen (Fig. 2) out of the eight resulted ones. For the authentication of the data and running a perfect supramolecular docking, the third binding pose was further subjected to Caver 3, to carry out quantum tunneling of the best active site binding pose [44]. In total, 24 quantum tunnels were developed by the Caver 3 from which the six most viable were chosen from the binding active site of the protein (Fig. 3A-F). The graph (Fig.  3G) depicts the quantitative analysis where the radius and tunnel length are dependent and independent variables, respectively, the radius aid in deducing the ligand's length, width, and height [45]. The tunnel was observed to be surrounded by amino acid residues, while the ligand paved its way through the tunnel (Fig. 3A). This very tunnel has been demonstrated in a graph with a radius of 2.0 Ǻ. The same phenomena were witnessed in Figure 3B and C, but with different radii of 1.5 Ǻ and 2.4 Ǻ, respectively [46]. Figure 3D shows a similar frame to 3B, where the ligand is in the tunnel center exhibiting a radius of 1.25 Å. Figure  3E shows a slightly different picture where the scattered amino acids are far away from the ligand. However, the radius' value of Figure 3E is identical to 3D, which is 1.25 Å. Lastly, in Figure 3F, fewer amino acids around the ligand depicted a radius with a value of 1.5 Å. By studying the quantitative analysis of the graph, it was concluded that Figure 3C tends to have the highest value which is 2.4 Å, and the lowest value is 1.25 Å found in Figure 3E and F.

Molecular optimization and docking
Molecular docking is a computational modeling technique that visualizes predictive pharmacophore complexes that occur between the ligand and the receptor protein via running them on program PyRx 0.8, a virtual screening tool which is based on Autodock Vina for the comprehension of the binding affinity between the ligand and the macromolecule when it is subjected to any condition depending on binding scores [47]. Before carrying out the molecular docking, UCSF Chimera Software Package (Version 1.14) was used to optimize the protein of interest and the five ligands. This software was also used to visualize the manner in which the binding positions if and when supramolecular docking will take place between the macromolecule and the ligand. Figure 4 shows the simultaneous qualitative (Fig. 4A) and quantitative (Fig. 4AB) forms of the control ligand serpentine complexed with the macromolecule hACE2 receptor. The qualitative form was developed via the UCSF Chimera Software, and here the predictive location of the hydrogen bonds between the specific amino acids are indicated using the magenta line. This was further conducted by the java interface-run program, LigPlot+ V.2.2 tool, to identify and visualize the hydrogen bonds and hydrophobic interactions between the ligand and the peripheral amino acid residues [32].
The protein was modeled in supramolecular docking with four ligands using the UCSF Chimera. The post-docking prediction of the binding poses on the ACE2 receptor through which the corresponding ligands bound to its protein, depending on the quantum tunneling mechanism conducted. The intermolecular spaces (red lines) and projected hydrogen bonds (magenta lines) are shown in Fig. 5.
The PyMOL-generated results (PDB format) were then subjected to visualization in the Ligplot+, which then displayed the number of hydrogen bond interactions and hydrophobic interactions of the ligand with the amino acids residing within its periphery (Fig. 6). Isovitexin (Fig.  6C) and apigenin (Fig. 6B) possessed the highest number of hydrogen bond interactions, while in the case of hydrophobic interactions, isovitexin tops apigenin. On the contrary, piperolactam A displayed the least number of hydrogen bonds along with hydrophobic interactions (Fig. 6A).
It was observed that the isovitexin complexed with the receptor possessed the highest binding affinity (-8.9 kcal/ mol) compared to the control ligand serpentine (−8.6 Kcal/ mol). This may be because the number of hydrogen bonds of the "isovitexin-1R4L" complex was higher than the "serpentine-1R4L" complex, whereas "apigenin-ACE2" protein had the least binding affinity value but was higher than the docked complex containing control ligand hydrogen bond, which plays a crucial role in protein-ligand binding by stabilizing the docked complex [48]. Besides hydrogen bonds, the higher number of hydrophobic interactions within the docked complex may also increase the binding affinity of the ligand toward the target protein [49,50]. Among the four test and control ligands, isovitexin has a higher number of carbonyl and hydroxyl moieties that may contribute the highest binding affinity of this ligand toward the protein.

Molecular dynamic simulation (20 ns)
In this research, Desmond (Desmond, Schrödinger, LLC, NY) has been used for molecular dynamic simulation (MDS) that was operated to investigate the receptor RMSD, receptor RMSF, Intra H-bonds, Cα, SASA, MolSA, Rg, and PSA for 20 ns. The RMSD values are considered to measure the average shift in a set of atoms for a specific frame with a reference frame [51]. This parameter governs a significant role in assisting the comparisons among different molecular structures and narrows down the extensive list of predictive conformations to a smaller set [52]. The greater the RMSD value, the less stable the docked complex during the simulation period, and vice versa [53]. In terms of the lowest RMSD values, all the experimental and control ligands demonstrated equal stability in their docked complexes. However, the stability pattern of the ligands, based on the RMSD values was piperolactam A > apigenin > serpentine > quercetin > isovitexin (Fig. 7A). From this trend, it was observed that the "piperolactam A-1R4L" complex showed the highest stability, whereas the "isovitexin-1R4L" complex exhibited the lowest stability during MDS analysis. The RMSF values were estimated for 20 ns to analyze the effect of a test ligand on 1R4L protein. RMSF value implies the denaturation tendency of a receptor protein at each point of the temporal trajectory (Fig. 7B). The higher the RMSF value, the lower the stability of the protein-ligand complex during the MDS and vice versa [53]. This phenomenon occurs due to a change in protein structure induced by ligand interaction when higher RMSF is attributed to the presence of tightly bonded structures like an alpha helix and beta-strands. At the same time, a lower RMSF refers to lose structures like coils, bends, and turns [54]. The ligand isovitexin exhibited the highest RMSF values than the control serpentine and the other experimental ligands, indicating lower stability due to enhanced flexibility. The control showed the least RMSF values to form the most stable complex. On the other hand, the docked complexes containing quercetin, piperolactam A, and apigenin showed lower stability than the control (Fig. 7B).
The hydrogen bond networks play a crucial role in strengthening the binding affinity of the protein and the ligand [55]. In our study, the "isovitexin-1R4L" complex showed the highest number of intramolecular hydrogen bond interactions among all ligands, including the control, indicating the most stable complex during the 20 ns of the simulation period. On the other hand, quercetin and apigenin exhibited a bit higher stability compared to control. The simulation findings demonstrated that intramolecular hydrogen bonds were absent in the docked complexes containing serpentine and piperolactam A, indicating the least stability in those complexes (Fig. 7C). The alpha-carbon atom is one of the mother parameters of MDS analysis which plays a vital role to get information on the motions of the protein-ligand complexes during simulation periods (Fig. 7D). Fastest (i.e., narrow region of alpha carbon atom) motions indicate the most biologically unstable conformations [56]. Therefore, the findings revealed that narrow regions of Cα were exhibited by the docked complex containing serpentine, apigenin, and piperolactam-A, indicating the fastest motion of their complexes, which resulted in the unstable conformations. However, the "isovitexin-1R4L" complex showed the widest region of carbon alpha, indicating stable conformations. Furthermore, the docked complex containing quercetin also exhibited a broader region of Cα compared to control, indicating the more stable conformation formed by the "quercetin-1R4L" complex than the serpentine-1R4L complex (Fig. 7D). There were no significant differences observed among serpentine, apigenin, and piperolactam A.
To better understand the effective interaction between the macromolecule and ligand, the SASA values are used to interpret the interactions between the surface of the docked complex and water molecules in which the proteinligand was submerged. For the stability of the proteins, hydrophobic interactions form between the non-polar amino acids through defending by polymer shielding in hydrophobic, which reduced SASA values [25,57]. The MDS analysis (Fig. 7E) revealed that serpentine (control) exhibited the most stability following the lowest SASA value. On the other hand, quercetin, and piperolactam A followed higher values of SASA, but isovitexin demonstrated the highest SASA value. The stability of a protein is also connected to its MolSA; consequently, a significant change in the MolSA of a complex might result in instability, which is highly undesirable [58,59]. After 20 ns of the simulation period, it was found that isovitexin showed the highest MolSA value among all five ligands (Fig. 7F). Thus, isovitexin can be considered the least stable and unfavorable in comparison to the other ligands. On the contrary, the most stable complex was formed between piperolactam A and 1R4L. Additionally, the docked complexes containing apigenin and quercetin were also displayed higher stability than the serpentine-1R4L complex (Fig. 7F).
The radius of gyration (Rg) is the benchmark to determine whether a structure has a stable, compact and folded conformation. The higher the Rg value, the more likely the ligand is flexible, thus possessing an unstable conformation. On the contrary, lower Rg values suggest a dense and closely packed structure [60]. Isovitexin exhibited the highest gyration radius, indicating the most flexible and unstable compound compared to the rest of the three experimental and control ligands. On the other hand, quercetin and apigenin exhibited higher stability than serpentine. From the simulation analysis, piperolactam-A was found to be the model ligand in terms of Rg (nm) value (Fig. 7G).
PSA is a crucial factor along with the lipophilicity to determine the ability of drugs to cross the blood-brain barrier [61]. The normal range of PSA for the substance (X) is 40 Å 2 < X ≤ 90 Å 2 , which ensures the best efficacy to permeate the blood-brain barrier, whereas out of this range is considered undesirable [62]. Only the control ligand can cross the blood-brain barrier among the five ligands since the PSA value of serpentine was to be found within the ideal range (Fig. 7H).
Although the research gained significant results in all aspects of "in silico molecular drug designing" for receptor specificity, few limitations also have been experienced. Firstly, certification of the ligands used in this study as targeted therapeutics cannot be carried out as long as no in vivo and clinical trial occurs. Secondly, the ligands of interest may show some mild to moderate side effects among the experimental subjects (animals like mice, rats, etc.). For instance, overuse of apigenin can cause stomach discomfort [63]; isovitexin can cause anxiety, affecting the brain cells [64]; serpentine shows hepatotoxicity sometimes [http://biosig.unimelb.edu.au/pkcsm/]; irritation and burning of the eyes, nose, throat, and skin in humans can emerge for piperolactam A; nausea, headache, and shortness of breath may sometimes result from quercetin overuse [https://www.rxlist.com/consumer_quercetin/ drugs-condition.htm]. Proper consideration of the effects of the dose is essential for these kinds of ligands.

Conclusion
This in silico study was conducted through a series of methodical procedures, which included: establishment of a phytochemical repository of 120 ligands based on a comprehensive literature review; selection of five top ligands including control based on drug likeliness through ADMET and QSAR-based ligand screening; in-depth analysis of the level of interaction among the receptor protein ACE2 (1R4L) and the respective test ligands through molecular docking of these ligands against the target receptor, optimized to a high-resolution crystalized level; and finally molecular dynamic simulation to attain the quantitative values for RMSD, RMSF, α-carbon, intra hydrogen bond, radius of gyration, SASA, MolSA, and PSA to observe the extent of the mobility of protein induced when in association with prospective drugs. Molecular docking analysis found all the test ligands to exhibit substantial binding affinity with the targeted macromolecule, while the Ligplot+ visualization demonstrated potential hydrogen and non-covalent bonds responsible for the strength of protein-ligand affinity. Moreover, the MDS suggested the hierarchy of five test ligands concerning least mobility of protein backbone, least fluctuation of protein alpha carbons, least flexibility of protein in complex with ligand, least available atomic van der Waals surface area, and protein surface area exposed to water for the protein-ligand complex, with a surface associated with heteroatoms and hydrogen atoms under desirable range. This study prognosticated apigenin, piperolactam A, and quercetin to be more stable and highly interactive flavonoid compounds. Thus, further in vivo investigation is recommended to determine their therapeutic potentiality and target specificity toward the target ACE2 receptor so that the SARS-CoV-2 influx can be prevented and the virus-derived heart injury can be protected. Although isovitexin exhibited the highest binding tendency with the hACE2 receptor, it was not suggested as the superior drug-like compound due to its low pharmacodynamics stability observed in MDS than the others.