Identification of Natural Lead Compounds against Hemagglutinin-Esterase Surface Glycoprotein in Human Coronaviruses Investigated via MD Simulation, Principal Component Analysis, Cross-Correlation, H-Bond Plot and MMGBSA

The pandemic outbreak of human coronavirus is a global health concern that affects people of all ages and genders, but there is currently still no effective, approved and potential drug against human coronavirus, as many other coronavirus vaccines have serious side effects while the development of small antiviral inhibitors has gained tremendous attention. For this research, HE was used as a therapeutic target, as the spike protein displays a high binding affinity for both host ACE2 and viral HE glycoprotein. Molecular docking, pharmacophore modelling and virtual screening of 38,000 natural compounds were employed to find out the best natural inhibitor against human coronaviruses with more efficiency and fewer side effects and further evaluated via MD simulation, PCA, DCCR and MMGBSA. The lead compound ‘Calceolarioside B’ was identified on the basis of pharmacophoric features which depict favorable binding (ΔGbind −37.6799 kcal/mol) with the HE(5N11) receptor that describes positive correlation movements in active site residues with better stability, a robust H-bond network, compactness and reliable ADMET properties. The Fraxinus sieboldiana Blume plant containing the Calceolarioside B compound could be used as a potential inhibitor that shows a higher efficacy and potency with fewer side effects. This research work will aid investigators in the testing and identification of chemicals that are effective and useful against human coronavirus.


Introduction
Human coronavirus (COVID-19) is a positive, single-stranded RNA virus that causes severe acute respiratory syndrome. It originated in China in 2019 and has spread to more than 210 countries [1]. The coronavirus outbreak appeared to be extremely dangerous

Structure Retrieval, Refinement and Evaluation
The X-ray crystallographic structure of Hemagglutinin-esterase (HE) surface glycoprotein with PDB ID 5N11 was retrieved from Protein Data Bank (PDB) (https://www.rcsb.org/ (accessed on 06 March 2022)), a freely available online database that contains the threedimensional structural data of macromolecules. The structure of the target protein was refined via BIOVIA Discovery studio [30] and evaluated via PROCHECK, which provides information about the stereochemistry of the protein structure via a Ramachandran Plot that describes the quality of the protein [31].

Selection of Ligands and Pharmacophore Generation
Fifteen antiviral ligands were retrieved against the target protein HE surface glycoprotein from the publicly accessible PubChem (https://pubchem.ncbi.nlm.nih.gov/ (accessed on 23 March 2022)) database [32] to generate a pharmacophore query. The use of pharmacophore is effective in computer-aided drug design (CADD). The pharmacophore model is an accumulation of common steric and electronic features that quickly filter through a huge number of a compound's library for a specific target to initiate or inhibit its biological response. The selected compounds were aligned and analyzed in terms of their chemical characteristics and common features were observed among them, such as hydrogen bond donor, acceptor, cationic, anionic, aromatic and hydrophobic [33]. Protein-ligand interactions were interpreted to achieve steric features. The pharmacophore of 15 active antiviral inhibitors was generated on the basis of RMSD values, common steric and chemical features, and a high binding affinity against Hemagglutinin esterase glycoprotein of human coronavirus using Ligand Scout 4.1.5. The Ligand Scout software rapidly generated a 3D pharmacophore from the structural data of small molecules in a fully automated and appropriate way [34].

Library Preparation and Virtual Screening
About 20 libraries that contain 38,000 natural compounds (Flavonoids, Traditional Chinese Medicinal compounds, highly selective inhibitors, antiviral and bioactive compounds) were downloaded from the SelleckChem (https://www.selleckchem.com/screening-libraries.html (accessed on 2 April 2022)) database on the basis of the Rule of Five and minimizing their energy via UCSF Chimera 1.14 before docking to make sure they were in the right conformation so the docking results would be more realistic. The system was subjected to energy minimization to start the production runs. The system was minimized using the steepest decent [35] and conjugate algorithm [36]. A total of 1500 steps of conjugate gradient algorithm were applied on the system in which after every 50th step, the deepest descent algorithm was applied. The freely available UCSF Chimera (https://www.cgl.ucsf.edu/chimera/ (accessed on 8 April 2022)) was used for visualization and analysis of molecular structures together with density maps, trajectories and sequence alignments. The pharmacophore model was used for the virtual screening of 38,000 natural compounds via Ligand Scout and it selected compounds that had the best pharmacophoric features and hit scores.

Docking Calculation and Interaction
After virtual screening, the top-20 compounds were selected on the basis of their pharmacophore fit score. Molecular docking of these compounds was performed through PyRx which contained Open Bebel, Vina Wizard, Autodock vina and python interpreter so that it could automatically convert files into the required file format. Active site residues were predicted via cocrystal structure and the CASTp server. The interpretation of H-bonds, polar, pi-anion, pi-alkyl, pi-donor hydrogen bond and hydrophobic interactions among HE receptor and studied compounds were visualized through UCSF Chimera 1.14 and PyMOL. On the basis of visualization, a best-hit compound was selected.

Lead Identification
Lead identification/optimization is an imperative step in drug design. All calculations (docking scores, interaction and ADMET analysis including MW, HBD, HBA, logP, PSA, rotatable bonds and rings) were achieved for the identification of a lead compound having the most suitable results by following the rules i.e., Ghose, Veber, Egan and rule of five (ROF). The compound with good interaction, best fit-score and binding affinity was selected as a potential inhibitor against HE surface glycoprotein.

Molecular Dynamic (MD) Simulations
The lead compound protein and reference complex were subjected to MD simulation to evaluate changes in the internal dynamics of the target protein. Amber tools were used to prepare input files while NAMD3 was utilized to conduct 100ns MD simulation with ff14SB and Gaff forcefields for protein (HE glycoprotein) and ligands (calceolarioside B and control), respectively. Both ligands and proteins were prepared via the Antechamber and Leap program of Amber tools while long-range electrostatic interactions were computed by the particle mesh Ewald method and short-range interactions, such as columbic and van der Waals interactions, were calculated with a cutoff of 10 Å. A specific number of Na + and Cl − counter ions and a TIP3P water box [41] of size 10 Å were introduced to imitate physiological salt concentration and to neutralize the whole system. A shake algorithm [42] was employed to constrain all bond lengths containing hydrogen bonds to heavy atoms while the particle mesh Ewald method [43] was utilized to calculate the longrange electrostatic interactions. The system was minimized using 10,000 steps and water equilibration was performed by using 10,000 steps. The temperature equilibrations were performed gradually at 200, 250, and 300 K temperatures for 5000 steps. After equilibration, the system was ready and prepared complexes were used to run MD simulation at a constant temperature of 310 K and 1 atm pressure. MD trajectories for both systems were analyzed to obtain RMSD, Rg, RMSF, energy decomposition, PCA, cross correlation and H-bond plot analysis.

Molecular Mechanics/Generalized Born Surface Area (MMGBSA) Analysis
Molecular Mechanics/Generalized Born Surface Area (MMGBSA) is the efficient force field technique to access binding free energy of a system (ligand-receptor) in kcal/mol [44]. Calculations were based on the last 300 frames and determined by using the following equations.
∆G bind = G complex − G protein − G ligand (1) ∆G gas = Bond + Angle + Dihed + EEL + VDWAAL (3) ∆G bind is the total binding free energy of the system (Equation (1)) which is calculated by employing Equation (2). T∆S is the change of conformational entropy on ligand binding at a given temperature. ∆G gas is the total of bond, angle, dihydral, EEL (electrostatic component of the internal energy) and van der Waals energy (Equation (3)). Internal energy is associated with the vibration and rotation of single bond torsional angles. Solvation free energy (∆G sol ) is the combination of ∆EGB (polar component of the solvation energy) and ∆ESURF (no polar component of the solvation energy) (Equation (4)).

Results
The target protein is the viral envelope protein, whose structure was retrieved from the Protein Data Bank with PDB ID 5N11 (Hemagglutinin esterase), 2.45 Å resolution and 0.249 Å R-free value. The HE receptor has two chains, 423 residues and 47,482 Da MW. The predicted structure was refined by removing the small compounds and water molecules via BIOVIA Discovery Studio, optimized, minimized and shown in Figure 1 along with the associated Ramachandran plot that provides information about the stereochemistry of the target protein. According to the graph, 87.2% residues are in most-favored regions, 12.1% are in additional allowed regions, while 0.7% are in generously allowed regions. Most of the residues are in most-favored regions and the overall quality of the HE protein is 95%, which indicates good quality structure.
ΔGgas = Bond + Angle + Dihed + EEL + VDWAAL (3) ΔGbind is the total binding free energy of the system (Equation (1)) which is calculated by employing Equation (2). TΔS is the change of conformational entropy on ligand binding at a given temperature. ΔGgas is the total of bond, angle, dihydral, EEL (electrostatic component of the internal energy) and van der Waals energy (Equation (3)). Internal energy is associated with the vibration and rotation of single bond torsional angles. Solvation free energy (ΔGsol) is the combination of ΔEGB (polar component of the solvation energy) and ΔESURF (no polar component of the solvation energy) (Equation (4)).

Results
The target protein is the viral envelope protein, whose structure was retrieved from the Protein Data Bank with PDB ID 5N11 (Hemagglutinin esterase), 2.45 Å resolution and 0.249 Å R-free value. The HE receptor has two chains, 423 residues and 47,482 Da MW. The predicted structure was refined by removing the small compounds and water molecules via BIOVIA Discovery Studio, optimized, minimized and shown in Figure 1 along with the associated Ramachandran plot that provides information about the stereochemistry of the target protein. According to the graph, 87.2% residues are in most-favored regions, 12.1% are in additional allowed regions, while 0.7% are in generously allowed regions. Most of the residues are in most-favored regions and the overall quality of the HE protein is 95%, which indicates good quality structure. According to the domain architecture of the HE protein, there are two domains, shown in Figure S1. The one is hema_esterase  while the second one is Hema HEFG (129-262), which is required for infection by recognizing the host cell receptor and helping with the fusion of the viral and host cell membrane.

Ligand-Based Virtual Screening and Molecular Docking
Fifteen antiviral compounds were used to generate the pharmacophore model which recognized the defined binding mode. The structures of these ligands along with their name and pharmacophore fit score are presented in Table 1. A total of 15 active According to the domain architecture of the HE protein, there are two domains, shown in Figure S1. The first one is hema_esterase  while the second one is Hema HEFG (129-262), which is required for infection by recognizing the host cell receptor and helping with the fusion of the viral and host cell membrane.

Ligand-Based Virtual Screening and Molecular Docking
Fifteen antiviral compounds were used to generate the pharmacophore model which recognized the defined binding mode. The structures of these ligands along with their name and pharmacophore fit score are presented in Table 1. A total of 15 active compounds were aligned via Ligand Scout and we generated a pharmacophore model by choosing the best features, such as HBD, HBA and aromaticity. Ligands with merged and selected pharmacophoric features are shown in Figure 2. The prepared pharmacophore model was used to screen a library of 38,000 natural compounds. Virtual screening has become a standard tool in drug discovery. After virtual screening the top-20, hits with the best pharmacophore fit score were selected for molecular docking with 5N11 receptor to explore their binding modes. Figure 3A,B illustrates docked poses within the active site and residues engaged in the binding. Docking energies of the top-20 ligands are shown in Table 2. The natural compound Calceolarioside B illustrates 11 interactions and was found most favorable for HE inhibition with the least binding energy of −7.8 kcal/mol. Polar amino acid residues i.e., Gln307, Cys311, Asn315, Asp289, Asp299 and nonpolar Phe313 created a hydrogen bond with OH of Calceolarioside B to inhibit the activity of the HE protein ( Figure 3C). The hydroxyl group (OH) increased the activity by participating in hydrogen bonding. The homovanillic-acid-HE complex formed two hydrogen bonds with OH of Phe313 while Gln307 formed hydrogen bond interactions with oxygen, respectively ( Figure 3D). As shown in Figure 3E, 2-(5-fluoro-2-methoxyphenyl)acetic acid formed a network of 02 conventional hydrogen bond interactions with polar Cys306 residue, and aromatic residue Phe313 while non-polar aliphatic Ala303, aromatic Tyr312 and polar uncharged Gln349 make carbon hydrogen bonds. ( Figure 3F) Hydroxytyrosol formed three conventional hydrogen bonds with Asn294, Trp292, and Arg296 while one carbon hydrogen bond with polar uncharged Ser298. Figure 3G depicts the hydrogen bond interactions of omarigliptin with Ser316, Arg291, Asn315, Ala303, Ser298, Asp299, Trp292, Asn293 and Gln307 residues. In addition, 4 -Methoxyresveratrol formed hydrogen bonds with Ala158, Tyr150, Tyr152 and Ala160 and pi-pi interactions with residue Ala173. The remaining 14 compounds displayed interactions with binding site residues as shown in Table S2. compounds were aligned via Ligand Scout and we generated a pharmacophore model by choosing the best features, such as HBD, HBA and aromaticity. Ligands with merged and selected pharmacophoric features are shown in Figure 2. The prepared pharmacophore model was used to screen a library of 38,000 natural compounds. Virtual screening has become a standard tool in drug discovery. After virtual screening the top-20, hits with the best pharmacophore fit score were selected for molecular docking with 5N11 receptor to explore their binding modes. Figure 3A,B illustrates docked poses within the active site and residues engaged in the binding. Docking energies of the top-20 ligands are shown in Table 2. The natural compound Calceolarioside B illustrates 11 interactions and was found most favorable for HE inhibition with the least binding energy of −7.8 kcal/mol. Polar amino acid residues i.e., Gln307, Cys311, Asn315, Asp289, Asp299 and nonpolar Phe313 created a hydrogen bond with OH of Calceolarioside B to inhibit the activity of the HE protein ( Figure 3C). The hydroxyl group (OH) increased the activity by participating in hydrogen bonding. The homovanillic-acid-HE complex formed two hydrogen bonds with OH of Phe313 while Gln307 formed hydrogen bond interactions with oxygen, respectively ( Figure 3D). As shown in Figure 3E, 2-(5-fluoro-2-methoxyphenyl)acetic acid formed a network of 02 conventional hydrogen bond interactions with polar Cys306 residue, and aromatic residue Phe313 while non-polar aliphatic Ala303, aromatic Tyr312 and polar uncharged Gln349 make carbon hydrogen bonds. ( Figure 3F) Hydroxytyrosol formed three conventional hydrogen bonds with Asn294, Trp292, and Arg296 while one carbon hydrogen bond with polar uncharged Ser298. Figure 3G depicts the hydrogen bond interactions of omarigliptin with Ser316, Arg291, Asn315, Ala303, Ser298, Asp299, Trp292, Asn293 and Gln307 residues. In addition, 4′-Methoxyresveratrol formed hydrogen bonds with Ala158, Tyr150, Tyr152 and Ala160 and pi-pi interactions with residue Ala173. The remaining 14 compounds displayed interactions with binding site residues as shown in Table S2. compounds were aligned via Ligand Scout and we generated a pharmacophore model by choosing the best features, such as HBD, HBA and aromaticity. Ligands with merged and selected pharmacophoric features are shown in Figure 2. The prepared pharmacophore model was used to screen a library of 38,000 natural compounds. Virtual screening has become a standard tool in drug discovery. After virtual screening the top-20, hits with the best pharmacophore fit score were selected for molecular docking with 5N11 receptor to explore their binding modes. Figure 3A,B illustrates docked poses within the active site and residues engaged in the binding. Docking energies of the top-20 ligands are shown in Table 2. The natural compound Calceolarioside B illustrates 11 interactions and was found most favorable for HE inhibition with the least binding energy of −7.8 kcal/mol. Polar amino acid residues i.e., Gln307, Cys311, Asn315, Asp289, Asp299 and nonpolar Phe313 created a hydrogen bond with OH of Calceolarioside B to inhibit the activity of the HE protein ( Figure 3C). The hydroxyl group (OH) increased the activity by participating in hydrogen bonding. The homovanillic-acid-HE complex formed two hydrogen bonds with OH of Phe313 while Gln307 formed hydrogen bond interactions with oxygen, respectively ( Figure 3D). As shown in Figure 3E, 2-(5-fluoro-2-methoxyphenyl)acetic acid formed a network of 02 conventional hydrogen bond interactions with polar Cys306 residue, and aromatic residue Phe313 while non-polar aliphatic Ala303, aromatic Tyr312 and polar uncharged Gln349 make carbon hydrogen bonds. ( Figure 3F) Hydroxytyrosol formed three conventional hydrogen bonds with Asn294, Trp292, and Arg296 while one carbon hydrogen bond with polar uncharged Ser298. Figure 3G depicts the hydrogen bond interactions of omarigliptin with Ser316, Arg291, Asn315, Ala303, Ser298, Asp299, Trp292, Asn293 and Gln307 residues. In addition, 4′-Methoxyresveratrol formed hydrogen bonds with Ala158, Tyr150, Tyr152 and Ala160 and pi-pi interactions with residue Ala173. The remaining 14 compounds displayed interactions with binding site residues as shown in Table S2. compounds were aligned via Ligand Scout and we generated a pharmacophore model by choosing the best features, such as HBD, HBA and aromaticity. Ligands with merged and selected pharmacophoric features are shown in Figure 2. The prepared pharmacophore model was used to screen a library of 38,000 natural compounds. Virtual screening has become a standard tool in drug discovery. After virtual screening the top-20, hits with the best pharmacophore fit score were selected for molecular docking with 5N11 receptor to explore their binding modes. Figure 3A,B illustrates docked poses within the active site and residues engaged in the binding. Docking energies of the top-20 ligands are shown in Table 2. The natural compound Calceolarioside B illustrates 11 interactions and was found most favorable for HE inhibition with the least binding energy of −7.8 kcal/mol. Polar amino acid residues i.e., Gln307, Cys311, Asn315, Asp289, Asp299 and nonpolar Phe313 created a hydrogen bond with OH of Calceolarioside B to inhibit the activity of the HE protein ( Figure 3C). The hydroxyl group (OH) increased the activity by participating in hydrogen bonding. The homovanillic-acid-HE complex formed two hydrogen bonds with OH of Phe313 while Gln307 formed hydrogen bond interactions with oxygen, respectively ( Figure 3D). As shown in Figure 3E, 2-(5-fluoro-2-methoxyphenyl)acetic acid formed a network of 02 conventional hydrogen bond interactions with polar Cys306 residue, and aromatic residue Phe313 while non-polar aliphatic Ala303, aromatic Tyr312 and polar uncharged Gln349 make carbon hydrogen bonds. ( Figure 3F) Hydroxytyrosol formed three conventional hydrogen bonds with Asn294, Trp292, and Arg296 while one carbon hydrogen bond with polar uncharged Ser298. Figure 3G depicts the hydrogen bond interactions of omarigliptin with Ser316, Arg291, Asn315, Ala303, Ser298, Asp299, Trp292, Asn293 and Gln307 residues. In addition, 4′-Methoxyresveratrol formed hydrogen bonds with Ala158, Tyr150, Tyr152 and Ala160 and pi-pi interactions with residue Ala173. The remaining 14 compounds displayed interactions with binding site residues as shown in Table S2. compounds were aligned via Ligand Scout and we generated a pharmacophore model by choosing the best features, such as HBD, HBA and aromaticity. Ligands with merged and selected pharmacophoric features are shown in Figure 2. The prepared pharmacophore model was used to screen a library of 38,000 natural compounds. Virtual screening has become a standard tool in drug discovery. After virtual screening the top-20, hits with the best pharmacophore fit score were selected for molecular docking with 5N11 receptor to explore their binding modes. Figure 3A,B illustrates docked poses within the active site and residues engaged in the binding. Docking energies of the top-20 ligands are shown in Table 2. The natural compound Calceolarioside B illustrates 11 interactions and was found most favorable for HE inhibition with the least binding energy of −7.8 kcal/mol. Polar amino acid residues i.e., Gln307, Cys311, Asn315, Asp289, Asp299 and nonpolar Phe313 created a hydrogen bond with OH of Calceolarioside B to inhibit the activity of the HE protein ( Figure 3C). The hydroxyl group (OH) increased the activity by participating in hydrogen bonding. The homovanillic-acid-HE complex formed two hydrogen bonds with OH of Phe313 while Gln307 formed hydrogen bond interactions with oxygen, respectively ( Figure 3D). As shown in Figure 3E, 2-(5-fluoro-2-methoxyphenyl)acetic acid formed a network of 02 conventional hydrogen bond interactions with polar Cys306 residue, and aromatic residue Phe313 while non-polar aliphatic Ala303, aromatic Tyr312 and polar uncharged Gln349 make carbon hydrogen bonds. ( Figure 3F) Hydroxytyrosol formed three conventional hydrogen bonds with Asn294, Trp292, and Arg296 while one carbon hydrogen bond with polar uncharged Ser298. Figure 3G depicts the hydrogen bond interactions of omarigliptin with Ser316, Arg291, Asn315, Ala303, Ser298, Asp299, Trp292, Asn293 and Gln307 residues. In addition, 4′-Methoxyresveratrol formed hydrogen bonds with Ala158, Tyr150, Tyr152 and Ala160 and pi-pi interactions with residue Ala173. The remaining 14 compounds displayed interactions with binding site residues as shown in Table S2. compounds were aligned via Ligand Scout and we generated a pharmacophore model by choosing the best features, such as HBD, HBA and aromaticity. Ligands with merged and selected pharmacophoric features are shown in Figure 2. The prepared pharmacophore model was used to screen a library of 38,000 natural compounds. Virtual screening has become a standard tool in drug discovery. After virtual screening the top-20, hits with the best pharmacophore fit score were selected for molecular docking with 5N11 receptor to explore their binding modes. Figure 3A,B illustrates docked poses within the active site and residues engaged in the binding. Docking energies of the top-20 ligands are shown in Table 2. The natural compound Calceolarioside B illustrates 11 interactions and was found most favorable for HE inhibition with the least binding energy of −7.8 kcal/mol. Polar amino acid residues i.e., Gln307, Cys311, Asn315, Asp289, Asp299 and nonpolar Phe313 created a hydrogen bond with OH of Calceolarioside B to inhibit the activity of the HE protein ( Figure 3C). The hydroxyl group (OH) increased the activity by participating in hydrogen bonding. The homovanillic-acid-HE complex formed two hydrogen bonds with OH of Phe313 while Gln307 formed hydrogen bond interactions with oxygen, respectively ( Figure 3D). As shown in Figure 3E, 2-(5-fluoro-2-methoxyphenyl)acetic acid formed a network of 02 conventional hydrogen bond interactions with polar Cys306 residue, and aromatic residue Phe313 while non-polar aliphatic Ala303, aromatic Tyr312 and polar uncharged Gln349 make carbon hydrogen bonds. ( Figure 3F) Hydroxytyrosol formed three conventional hydrogen bonds with Asn294, Trp292, and Arg296 while one carbon hydrogen bond with polar uncharged Ser298. Figure 3G depicts the hydrogen bond interactions of omarigliptin with Ser316, Arg291, Asn315, Ala303, Ser298, Asp299, Trp292, Asn293 and Gln307 residues. In addition, 4′-Methoxyresveratrol formed hydrogen bonds with Ala158, Tyr150, Tyr152 and Ala160 and pi-pi interactions with residue Ala173. The remaining 14 compounds displayed interactions with binding site residues as shown in Table S2.          ADMET properties and heatmap for toxicity analysis were determined for the top-six compounds in Table 3 and Figure 4, respectively, that lay within the acceptable toxicity profile. The bioavailability radar of the six top hits is shown in Figure 5. Phytochemical calceolarioside B was selected as a lead compound against HE glycoprotein of human coro-navirus on the basis of its best binding affinity, RMSD, good interactions, pharmacophore fit score and ADMET properties.   ADMET properties and heatmap for toxicity analysis were determined for the topsix compounds in Table 3 and Figure 4, respectively, that lay within the acceptable toxicity profile. The bioavailability radar of the six top hits is shown in Figure 5. Phytochemical calceolarioside B was selected as a lead compound against HE glycoprotein of human coronavirus on the basis of its best binding affinity, RMSD, good interactions, pharmacophore fit score and ADMET properties.

Molecular Dynamics Simulations
To further analyze the stability of the lead compound-receptor and control complex, MD simulation was performed, and we evaluated MD trajectories and determined the RMSD, RMSF, Radius of Gyration, principal component analysis (PCA), cross correlation, no. of hydrogen bonds and MMGBSA.  Figure 6A). The RMSD value of the control and the complex did not vary significantly, staying nearly constant over the course of the simulation and represented rigidity ( Figure 6B). To identify the dynamic behavior of most mobile residues and the effect of calceolarioside B binding on the flexibility of the target protein, RMSF was investigated. According to the RMSF plot, protein residues did not experience

Molecular Dynamics Simulations
To further analyze the stability of the lead compound-receptor and control complex, MD simulation was performed, and we evaluated MD trajectories and determined the RMSD, RMSF, Radius of Gyration, principal component analysis (PCA), cross correlation, no. of hydrogen bonds and MMGBSA.  Figure 6A). The RMSD value of the control and the complex did not vary significantly, staying nearly constant over the course of the simulation and represented rigidity ( Figure 6B). To identify the dynamic behavior of most mobile residues and the effect of calceolarioside B binding on the flexibility of the target protein, RMSF was investigated. According to the RMSF plot, protein residues did not experience much flexibility when binding with a hit and reference compound but a noteworthy escalation in the flexibility of amino acids was analyzed between 165-205 residues while an escalation in flexibility was remarkable in the case of the control. Except for the specific region (165-205), both systems illustrated similarity in residual fluctuations as the average RMSF value of the complex was 1.10 ± 0.25 Å and 1.50 ± 0.30 Å for the control. The higher RMSF value at the C and N terminal and the middle residues represented loop regions that fluctuated more than other secondary structures while a lower RMSF represented relatively rigid residues. We detected the compactness of hemagglutinin esterase by approaching the radius of gyration (Rg) of carbon alpha atoms. It provided insights into the overall protein dimensions and enabled evaluations of the modifications to the tertiary structure of the protein throughout the simulation. Figure 6C displays a fluctuation near 40 and 58 ns in the complex system while the control represents significant fluctuations at 35, 68 and 98 ns with an abrupt decrease near 60 ns and a sustained average 20.50 Å Rg value. Overall, there was no major variation in the complex Rg throughout the simulation which showed that there were no unfolding events or loose packing, and revealed the extremely compacted nature of protein-ligand complex, and the complex maintained a 21.74 ± 0.18 Å Rg value.

System Stability, Fluctuation and Radius of Gyration
Biomedicines 2023, 11, x 11 of 17 escalation in the flexibility of amino acids was analyzed between 165-205 residues while an escalation in flexibility was remarkable in the case of the control. Except for the specific region (165-205), both systems illustrated similarity in residual fluctuations as the average RMSF value of the complex was 1.10 ± 0.25 Å and 1.50 ± 0.30 Å for the control. The higher RMSF value at the C and N terminal and the middle residues represented loop regions that fluctuated more than other secondary structures while a lower RMSF represented relatively rigid residues. We detected the compactness of hemagglutinin esterase by approaching the radius of gyration (Rg) of carbon alpha atoms. It provided insights into the overall protein dimensions and enabled evaluations of the modifications to the tertiary structure of the protein throughout the simulation. Figure 6C displays a fluctuation near 40 and 58 ns in the complex system while the control represents significant fluctuations at 35, 68 and 98 ns with an abrupt decrease near 60ns and a sustained average 20.50 Å Rg value. Overall, there was no major variation in the complex Rg throughout the simulation which showed that there were no unfolding events or loose packing, and revealed the extremely compacted nature of protein-ligand complex, and the complex maintained a 21.74 ± 0.18 Å Rg value.
( Figure 6D) The time evolution plot of H-bonds determines the formation and stability of hydrogen bonds throughout the simulation time as H-bonds play a significant role in drug specificity, metabolism and absorption [45,46]. The results illustrated that calceolarioside B formed up to four H-bonds with 88.34% occupancy and depicted the stable nature of the complex while the control formed up to ten hydrogen bonds with a minimum of six H-bonds throughout the simulation period.

Principal Component Analysis (PCA)
A Principal Component Analysis (PCA) was employed to detect the protein's conformational changes mediated by calceolarioside B binding and reveal the collective motions of MD trajectories. According to Figure 7, PC1, PC2, PC3 and eigenvalues of receptor ( Figure 6D) The time evolution plot of H-bonds determines the formation and stability of hydrogen bonds throughout the simulation time as H-bonds play a significant role in drug specificity, metabolism and absorption [45,46]. The results illustrated that calceolarioside B formed up to four H-bonds with 88.34% occupancy and depicted the stable nature of the complex while the control formed up to ten hydrogen bonds with a minimum of six H-bonds throughout the simulation period.

Principal Component Analysis (PCA)
A Principal Component Analysis (PCA) was employed to detect the protein's conformational changes mediated by calceolarioside B binding and reveal the collective motions of MD trajectories. According to Figure 7, PC1, PC2, PC3 and eigenvalues of receptor was plotted against the respective eigenvector index for the first 20 modes of motion. PC analysis indicated conformational changes in all clusters where the blue region exhibited the most significant movements, the white region represented intermediate movements and the red region displayed the least flexible movements. Overall protein movement was controlled by eigenvectors, especially the higher ones and the top-five eigenvectors in our system demonstrated dominant movements with eigenvalues of 18.0-59.9% while the remaining eigenvectors had lower eigenvalues. According to the PCA plot, the PC1 cluster retained the highest variability of 17.98%, PC2 illustrated 10.63% variability, while PC3 showed minimal variability (8.09%). The minimal variability of PC3 indicates highly stabilized protein ligand binding and a compact structure when compared to the PC1 and PC2 clusters.

Positive-Negative Correlation Movements of Residues
Dynamic cross-correlation maps represent inter residual motions computed via MD trajectories (Figure 8). The cyan and magenta color depicts strongly correlated (positive) and anticorrelated (negative) motions, respectively, between essential residues throughout the MD simulation. The correlated residues were more than 0.8 while anticorrelated residues were <−0.4. Positive correlation confirms RMSD and revealed a high stability. As shown in Figure 8A, the control (standard compound-HE) depicts a positive correlation and also depicts some of the negatively correlated movements ( Figure 8B). The lead compound and receptor are significantly correlated, and positively correlated movements were extremely notable at residues 100-130 and 250-350 (active site region) and a higher percentage of pairwise-correlated residues represent the stable binding of calceolarioside B with the HE protein.

Positive-Negative Correlation Movements of Residues
Dynamic cross-correlation maps represent inter residual motions computed via MD trajectories (Figure 8). The cyan and magenta color depicts strongly correlated (positive) and anticorrelated (negative) motions, respectively, between essential residues throughout the MD simulation. The correlated residues were more than 0.8 while anticorrelated residues were <−0.4. Positive correlation confirms RMSD and revealed a high stability. As shown in Figure 8A, the control (standard compound-HE) depicts a positive correlation and also depicts some of the negatively correlated movements ( Figure 8B). The lead compound and receptor are significantly correlated, and positively correlated movements were extremely notable at residues 100-130 and 250-350 (active site region) and a higher percentage of pairwise-correlated residues represent the stable binding of calceolarioside B with the HE protein.

Binding Energy Landscape and Energy Decomposition Analysis
To estimate the contribution of individual residues towards HE protein's inhibition, MMGBSA and energy decomposition analyses were performed. Per-residue energy decomposition analysis showed the contribution of different amino acids to the overall binding energy. According to the energy decomposition graph, the highest contributing residues Val304, Phe313, Gln349 and Asn300 with −3.3, −2.5, −2, −1.5 kcal/mol energies interacted with the ligand and are highlighted in Figure 9A.
Highly dynamic, cost-effective and computer-derived MMGBSA analysis computes the binding free energy of the protein-ligand complex at the molecular level that might be extremely beneficial for drug design ( Figure 9B). The binding free energy (ΔGbind) of calceolarioside B complex with HE protein is −37.6799 kcal/mol. Data reveal that van der Waals interactions (VDWAALS) significantly contribute (−46.4165 kcal/mol) to the binding free energies while EGB was 23.0200 kcal/mol and EEL was −9.4782 kcal/mol. The ΔGgas (bond + angle + dihed + EEL + VDWAALS) was the highest energy with −55.8948 kcal/mol value (Table S1). The Calceolarioside B-HE complex represented the lowest negative values, indicating stability and favorable binding of calceolarioside B in the active site of the HE receptor.

Discussion
The pandemic outbreak of novel human coronavirus spread into several other countries. In 2020, WHO declared a global health emergency based on the growing number of

Binding Energy Landscape and Energy Decomposition Analysis
To estimate the contribution of individual residues towards HE protein's inhibition, MMGBSA and energy decomposition analyses were performed. Per-residue energy decomposition analysis showed the contribution of different amino acids to the overall binding energy. According to the energy decomposition graph, the highest contributing residues Val304, Phe313, Gln349 and Asn300 with −3.3, −2.5, −2, −1.5 kcal/mol energies interacted with the ligand and are highlighted in Figure 9A.

Binding Energy Landscape and Energy Decomposition Analysis
To estimate the contribution of individual residues towards HE protein's inhibition, MMGBSA and energy decomposition analyses were performed. Per-residue energy decomposition analysis showed the contribution of different amino acids to the overall binding energy. According to the energy decomposition graph, the highest contributing residues Val304, Phe313, Gln349 and Asn300 with −3.3, −2.5, −2, −1.5 kcal/mol energies interacted with the ligand and are highlighted in Figure 9A.
Highly dynamic, cost-effective and computer-derived MMGBSA analysis computes the binding free energy of the protein-ligand complex at the molecular level that might be extremely beneficial for drug design ( Figure 9B). The binding free energy (ΔGbind) of calceolarioside B complex with HE protein is −37.6799 kcal/mol. Data reveal that van der Waals interactions (VDWAALS) significantly contribute (−46.4165 kcal/mol) to the binding free energies while EGB was 23.0200 kcal/mol and EEL was −9.4782 kcal/mol. The ΔGgas (bond + angle + dihed + EEL + VDWAALS) was the highest energy with −55.8948 kcal/mol value (Table S1). The Calceolarioside B-HE complex represented the lowest negative values, indicating stability and favorable binding of calceolarioside B in the active site of the HE receptor.

Discussion
The pandemic outbreak of novel human coronavirus spread into several other countries. In 2020, WHO declared a global health emergency based on the growing number of cases and a situation that was growing worse on a daily basis. To handle this situation, it Highly dynamic, cost-effective and computer-derived MMGBSA analysis computes the binding free energy of the protein-ligand complex at the molecular level that might be extremely beneficial for drug design ( Figure 9B). The binding free energy (∆G bind ) of calceolarioside B complex with HE protein is −37.6799 kcal/mol. Data reveal that van der Waals interactions (VDWAALS) significantly contribute (−46.4165 kcal/mol) to the binding free energies while EGB was 23.0200 kcal/mol and EEL was −9.4782 kcal/mol. The ∆G gas (bond + angle + dihed + EEL + VDWAALS) was the highest energy with −55.8948 kcal/mol value (Table S1). The Calceolarioside B-HE complex represented the lowest negative values, indicating stability and favorable binding of calceolarioside B in the active site of the HE receptor.

Discussion
The pandemic outbreak of novel human coronavirus spread into several other countries. In 2020, WHO declared a global health emergency based on the growing number of cases and a situation that was growing worse on a daily basis. To handle this situation, it is necessary to develop new drugs to treat COVID-19. Therefore, we used a multidisciplinary field, computer-aided drug design widely used to find new drug candidates in less time and at a reduced cost.
In this research, we made use of different bioinformatics tools to find a natural inhibitor against the HE surface glycoprotein of human coronavirus. In recent years, the use of natural compounds against viral infections has been found to be effective and is gaining importance. Natural compounds are less toxic, and less harmful to human health and are being analyzed to understand whether they could inhibit coronavirus. This virus shares sequence similarity with beta coronaviruses which possess the HE protein that interacts with various types of sialic acid, removes acetyl groups from O-acetylated sialic acid and play a role in binding to the target cell.
Ligands and receptor were obtained from PubChem and PDB respectively, then minimize their energy in order to reduce the overall potential of the receptor and ligands and to make sure that they were in the right conformation with low delta G values so as to be considered close to the biological system. Libraries of natural inhibitor compounds were downloaded from the SelleckChem database which provided the antiviral, antifungal and anti-inflammatory effects. A pharmacophore model was generated on the basis of shared steric and electronic features of all known active and antiviral compounds with a wide range of structural diversity and activities were aligned which were responsible for the biological interactions. The pharmacophore model explains how structurally distinct ligands can bind to the same side of the receptor.
Virtual screenings of thousands of natural compounds have been performed to discover novel molecules from a library of 38,000 compounds, docked via PyRx and to determine the interaction between the small molecule and the active site of the target protein at an atomic level. Resultantly, we selected the best compound, calceolarioside B, on the basis of its best docking affinity, RMSD and other physiochemical properties. The lead compound tightly bound with Asn315, Val304, Phe313, Asp289, Cys311, Asn300 and Gln307 amino acids and stabilized the active site of the HE receptor. Calceolarioside B is derived from the roots and leaves of Fraxinus sieboldiana Blume plant, which is member of the Oleaceae family, commonly known as the ash tree that is found in various regions of the world. It is native to China Southeast, Japan and Korea. In northern areas of Pakistan, Fraxinus sieboldiana Blume plant is usually used to treat malaria and pneumonia. Metabolites and extracts from this plant exhibit a wide range of biological actions, including anticancer, anti-inflammatory, antioxidant, antimicrobial, hepatoprotective, antiallergic and anti-viral properties.
A toxicity analysis evaluated the safety of the potential drug candidate, which indicated that calceolarioside B is safe to use in the future against the HE glycoprotein of human coronavirus. In checking for toxicity, the Lipinski rule was required to be followed to check whether the drug was toxic or non-toxic. The Calceolarioside B compound violated two Lipinski rules, but a drug or compound with two violations is acceptable while more than two violations are not acceptable. RMSD, which computes the average distance and the binding of calceolarioside B, revealed stability in the HE receptor. The RMSF graph did not represent major fluctuations in the target protein after binding of a hit compound. The radius of gyration value represents the compactness and stabilized folding in the phytochemical bound complex. Calceolarioside B formed a lot of interactions with the HE receptor, but the H-bond played a significant role by stabilizing the complex. PCA, DCCR and MMGBSA represent compactness and stability in the lead compound-HE receptor. Generally, the analysis of the MD simulation trajectory revealed the stable and energetically favorable complex formation in the presence of a lead phytochemical. These findings and their implications are discussed in the broadest context possible. This study will help researchers to evaluate the compounds that are effective and beneficial against human coronavirus.

Conclusions
HE glycoprotein (5N11) is involved in causing COVID-19 disease in humans. This research work was designed to find a natural compound that can act as inhibitor against HE glycoprotein of human coronavirus within a reduced time and cost by using different bioinformatics tools. Pharmacophore modeling, virtual screening and molecular docking helped to filter out calceolarioside B as having a low binding affinity with the target protein.
Fraxinus sieboldiana Blume, a medicinal plant with rich phytochemical compounds, out of which calceolarioside B is one of the compounds that show antiviral activity, inhibits the replication of coronavirus, stabilizing the structure and energy of the HE receptor indicated via MD simulation and MMGBSA analysis. The selected lead phytochemical must be validated in future through in vitro and in vivo studies. It is concluded that Calceolarioside B that is present in the root bark and leaves of the Fraxinus sieboldiana Blume plant, is an effective lead compound in the case of novel coronavirus.

Supplementary Materials:
The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/biomedicines11030793/s1. Figure S1: Domain architecture of HE Protein. Table S1. Calculated values of MMGBSA with different parameters of HE proteincalceolarioside B complex. Table S2 Funding: Article processing charges for this manuscript were supported by the University of Oradea, Romania.