Natural inhibitors for severe acute respiratory syndrome coronavirus 2 main protease from Moringa oleifera, Aloe vera, and Nyctanthes arbor-tristis: molecular docking and ab initio fragment molecular orbital calculations

The novel coronavirus 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has spread rapidly worldwide, and new drug treatments for COVID-19 are urgently required. To find the potential inhibitors against the main protease (Mpro) of SARS-CoV-2, we investigated the inhibitory potential of naturally occurring compounds from the plants Moringa oleifera, Aloe vera, and Nyctanthes arbor-tristis, using molecular docking, classical molecular mechanics optimizations, and ab initio fragment molecular orbital (FMO) calculations. Of the 35 compounds that we simulated, feralolide from Aloe vera exhibited the highest binding affinity against Mpro. Therefore, we proposed novel compounds based on the feralolide and investigated their binding properties to Mpro. The FMO results indicated that the introduction of a hydroxyl group into feralolide significantly enhances its binding affinity to Mpro. These results provide useful information for developing potent Mpro inhibitors. Supplementary Information The online version contains supplementary material available at 10.1007/s11224-022-02021-y.


Introduction
The novel coronavirus disease 2019 (COVID-19) responsible for the ongoing pandemic was first reported in December 2019 [1]. Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the causative agent of COVID-19, has spread rapidly around the world [2,3]. It is an enveloped positive-sense single-stranded RNA virus and belongs to the Coronaviridae family, which causes a wide range of infections [4]. It contains one of the coronavirus nonstructural proteins, which is known as the main protease (Mpro). Mpro cleaves the viral polyproteins, pp1a, and pp1ab, into 16 nonstructural proteins that are important for viral replication and maturation. Moreover, Mpro plays an important role in virus entry into host cells. Therefore, Mpro is a potential target for drugs to treat coronavirus infection because inhibition of Mpro halts virus entry and subsequent replication [5]. Currently, there is no effective antiviral therapy for COVID-19. Although vaccines have been developed against it, their safety and long-term effectiveness are still unclear [6], and no vaccine is 100% protective. To address these problems, we focused on identifying natural compounds that inhibit the Mpro of SARS-CoV-2 as potential drug leads based on the ab initio molecular simulations.
Medicinal plants are a valuable source of therapeutic agents and new leads for drug discovery [6,7]. Approximately, 85% of the world's population depends on medicinal plants for primary health care needs, and 80% of all synthetic drugs are derived from compounds found in medicinal plants [7]. Natural medicines are effective, safe, and affordable and have fewer side effects compared with synthetic drugs [8].
Aloe vera, also called Aloe barbadensis Miller, a species within the Asphodelaceae family, is a very short-stemmed, perennial, xerophytic, and succulent plant with triangular, green fleshy leaves [9,10]. It is the most common aloe variety and contains a wide diversity of bioactive compounds including terpenoids, lectins, flavonoids, anthraquinones, fatty acids, polysaccharides, sterols, tannins, enzymes, vitamins, and minerals [8,11]. It has previously been reported to have antiviral, antimicrobial, anti-inflammatory, and antidiabetic effects [8,12,13]. The antiviral activity of Aloe vera is due to its content of anthraquinones, flavonoids, minerals, vitamins, phenolic acids, sterols, and polysaccharides [8,13]. Aloe-emodin (A10 of Fig. 1), an anthraquinone present in Aloe vera, has been reported to have antiviral activity against the following viruses: herpes simplex types 1 and 2, varicella-zoster, pseudorabies, and influenza [13,14]. Substances extracted from Aloe vera leaves have antiviral effects against human cytomegalovirus in vitro [15]. In addition, Rezazadeh et al. have investigated the anti-herpes simplex virus type 1 activity of Aloe vera gel in in vitro cell culture [16]. Recently, Gansukh et al. have reported the ultrasound-mediated anti-influenza activity of the anthraquinones aloin and aloe-emodin (A10) from Aloe vera [17].
Nyctanthes arbor-tristis, commonly known as night jasmine, is a traditional medicinal plant that belongs to the family Oleaceae [18][19][20]. This plant is native to India and found in gardens as an ornamental plant with high medicinal value in Ayurveda, a form of alternative medicine originating in India. It is widely distributed in the southern sub-Himalayan regions, as far as Godavari [19,20]. Flavonoids, steroids, terpenes, alkaloids, and aliphatic compounds have been isolated from different tissues of Nyctanthes arbor-tristis. The main classes of medicinally interesting chemicals produced by this plant are alkaloids and glycosides [20]. Various studies on the pharmacological action of Nyctanthes arbortristis extracts have been reported, including antibacterial, anti-diabetic, antioxidant, anti-cancer, anti-allergic, antimalarial, and antiviral effects. Ethanolic extracts, n-butanol fractions, and two pure compounds, arbortristosides A and C, isolated from Nyctanthes arbor-tristis have shown inhibitory activity against encephalomyocarditis virus and Semliki Forest virus [20][21][22].
In previous molecular simulations [23] based on the molecular docking and ab initio fragment molecular orbital (FMO) methods, we investigated the binding properties of compounds extracted from Moringa oleifera to the Mpro of SARS-CoV-2. Moringa oleifera is a plant belonging to the Moringaceae family and has antifungal, antioxidant, antibacterial, anti-inflammatory, diuretic, and hepato-protective properties. Of the 12 natural products extracted from Moringa oleifera, the compound niaziminin (M9 shown in Fig. S1 of Supplementary material) has been found to have the highest binding affinity to the Mpro. Therefore, we have proposed novel compounds based on the structure of niaziminin and examined their binding properties to the Mpro, revealing that the addition of a hydroxyl group to niaziminin enhances its binding affinity to the Mpro.
The present study aims to extend this work to Aloe vera and Nyctanthes arbor-tristis to find novel compounds with higher binding affinities to the Mpro. To precisely investigate the binding properties of the natural products to the Mpro, we here performed a combination of approaches, including molecular docking, molecular mechanics (MM) optimization, and ab initio FMO calculations. The results revealed that feralolide (A14 of Fig. 1) from Aloe vera exhibits the highest binding affinity to the Mpro. Therefore, it is likely to be a potent Mpro inhibitor and a drug lead compound. Next, we proposed the novel compounds based on feralolide (A14) and investigated their binding properties. In addition, we conducted molecular dynamics (MD) simulations to confirm the stability of our proposed compounds at the ligand-binding pocket of the Mpro. Based on the FMO and MD results, we proposed some derivatives of natural compounds as candidates for potent inhibitors against the Mpro of SARS-CoV-2.
The structures of the selected compounds were fully optimized using the B3LYP/6-31G (d, p) method of the ab initio molecular orbital calculation program Gaussian 09 (G09) [34]. The charge distribution of each optimized structure was evaluated using the restrained electrostatic potential (RESP) analysis [35] using G09 [34] with the HF/6-31G (d) method, and the RESP charges were used as the charge parameters in the MM force fields of the compounds. These RESP charges are essential for the docking simulations, MM optimizations of the Mpro-ligand complexes, to precisely describe the electrostatic interactions between the compounds and the Mpro.
The 3D structure of the Mpro [36] was downloaded from the Protein Data Bank (PDB) as PDB ID: 6LU7 [37]. This protein has seven histidine (His) residues, each of which can exist in three different protonation states (i.e., Hid, Hie, and Hip) depending on the acidity of the surrounding environment. When the acidity of the environment is below pH6, a His residue tends to be in the protonated Hip state, whereas other His residues in the same protein could still be in the Hid or Hie state depending on the local structure around them. We used the PROPKA3.1 software [38,39] to predict the exact protonation state of the His residues of Mpro. Of the seven His residues, only His64 was assigned as Hip; His41, His80, His163, and His172 were assigned to the Hid state, whereas His164 and His246 were assigned to the Hie state considering the local structure around these His residues. The ionized state was adopted for the remaining ionizable amino acid residues in the Mpro.

Molecular docking of the compounds to the main protease and MM optimization
To find the best inhibitors in the Mpro binding site, we conducted molecular docking simulations using the program AutoDock4.2.6 [40]. In the docking simulations, the size of the grid box was set to 19.5 × 19.5 × 19.5 Å 3 , which was almost 1.5 times the size of the compound, and the center of the grid box was set to the center of the compound in the co-crystallized Mpro-compound complex (PDB ID: 6LU7). The number of created candidate poses was 256, and the threshold distance for clustering these poses was set as 1.5 Å. The maximum number of energy evaluations (ga_num_evals) for each run was set as 2,500,000. From the various clusters generated by AutoDock4.2.6 [40], we selected the three clusters with the largest number of poses, and the representative structures of these clusters were used in the subsequent MM and FMO calculations.
To obtain stable structures for the Mpro-ligand complexes, the representative structures of the clusters obtained from the docking simulations were fully optimized in water using the classical MM method. In the MM optimizations, approximately 1800 water molecules existing within 8 Å of the surface of the complex were explicitly considered. The MM and MD simulation program AMBER12 [41] was used. The AMBER FF99-SBLIN force field [42], TIP3P model [43], and general AMBER force field [44] were assigned to the Mpro, water molecules, and compounds, respectively. The criterion for the convergence of structure optimization was set as 0.0001 kcal/mol/Å.

Ab initio FMO calculations for the main proteaseligand complexes
To clarify the specific interactions and binding affinity of each compound to the Mpro, we investigated the electronic properties of the Mpro-ligand complex with explicit waters using the ab initio FMO method [45]. This method has obtained accurate results, comparable with those from experiments, for many biomolecules. As water molecules can make specific contributions to the interactions between the Mpro and a compound, water molecules existing within 10 Å of each compound were considered explicitly. The number of water molecules used in the FMO calculations was approximately 150 for all clusters. To predict the binding affinity of each compound to the Mpro, we evaluated the total inter-fragment interaction energies (IFIEs) [46] between the compounds and the Mpro residues using the ab initio FMO method [45].
In the present study, we did not consider the effect of entropy on the binding affinity because a vibrational analysis for the solvated Mpro-ligand complex would not have been practical when using the ab initio FMO method. In addition, the entropic effect was unlikely to be markedly different for each of the compounds, as they had similar chemical structures and would bind to the same site in the Mpro. Therefore, we investigated the total IFIEs between the Mpro residues and each compound using the ab initio FMO calculations and estimated the trend of binding affinity, assuming that the entropic effect was the same for each compound.
In the FMO calculations, the MP2/6-31G (d) method [47,48] of the ABINIT-MP version 6.0 FMO calculation program [49] was used. Each amino acid of the Mpro, compound, and each water molecule were assigned as fragments in the FMO calculations. This fragmentation enabled us to analyze the interactions between each Mpro residue and the compound, including any effects due to the solvating water molecules.
In our previous study [50], we investigated the binding of ligands to the androgen receptor protein using the same FMO calculations. The evaluated total IFIEs of the ligands with the androgen receptor residues were confirmed to correlate well with the binding affinities of these ligands obtained experimentally. The coefficient of determination was 0.94 using all nine of the ligands studied, confirming that our evaluated total IFIEs correlated well with the observed binding affinities of these ligands for the androgen receptor. Therefore, the present FMO calculations are expected to estimate the binding affinity of the compounds and the Mpro accurately.

Optimized structures of the main protease-ligand complexes
In our previous molecular simulations [23], we considered 12 natural compounds extracted from Moringa oleifera and investigated their binding affinities to the Mpro. Their chemical properties and structures are shown in Table S1 and Fig. S1 of the Supplementary material (SM). In addition, in the present study, we considered 14 compounds from Aloe vera and 9 compounds from Nyctanthes arbor-tristis as candidate inhibitors against Mpro. Their pharmacokinetic properties evaluated using the SwissADME web tool [33] are listed in Tables S2 and S3  At first, the structures of the compounds were optimized, and the RESP charges were evaluated by G09 [34]. Using the optimized structure and RESP charge, each compound was docked into the binding site of the Mpro using the AutoDock4.2.6 program [40]. The candidate structures of Table 1 This table shows the compounds extracted from Aloe vera with their lowest binding energy (BE) (kcal/mol), the number of poses, and the Mpro residues involved in H-bonds for the selected clusters obtained using the AutoDock4.2.6 program [40]. The 256 poses created were classified into several clusters based on the structural similarity, and each cluster was ranked in the order of the BE between the Mpro and the compound. We selected three clusters with the highest number of poses and evaluated the total inter-fragment interaction energy (IFIE) (kcal/ mol) of each compound with all Mpro residues using the FMO method. These values are listed in the last column the Mpro-ligand complex produced by the AutoDock4.2.6 program [40] were grouped into several clusters according to their structural similarities, and the clusters were ranked based on the binding energy between the Mpro and the compound. We selected three clusters with the largest number of poses because there was a higher probability that the compounds possessed one of the candidate structures in these clusters. Therefore, we chose representative structures from these clusters and optimized them with explicit water molecules using the MM method of AMBER12 [41]. For the optimized structures, the total IFIEs between the compounds and each of the Mpro residues were precisely evaluated using the ab initio FMO method [45] to determine which compound bound most strongly to the Mpro.
The results of the docking, MM, and FMO calculations for the compounds from Aloe vera and Nyctanthes arbortristis are listed in Tables 1 and 2, respectively, where the lowest binding energy (BE) in kcal/mol, number of candidate poses, and Mpro residues involved in hydrogen bonds with each compound are included for the selected clusters. In the last column of these Tables, the total IFIE in kcal/mol for each compound with all Mpro residues, evaluated using the ab initio FMO method is listed. Notably, Tables 1 and 2 indicate that the top-ranked cluster created by AutoDock4.2.6 did not necessarily correspond to the most stable structure in the FMO calculations. For example, as listed in Table 1, cluster 4 of compound A1 had the highest total IFIE, indicating that A1 was bound most strongly to the Mpro in the conformation of cluster 4, not in cluster 1. Therefore, Tables 1 and 2 suggested that the selection of clusters from the AutoDock4.2.6 results had to be done with caution. In the following analysis, we employed the cluster with the highest total IFIE.

Binding affinities of the compounds and the main protease
The total IFIEs between the Mpro residues and the natural compounds extracted from Moringa oleifera [23], Aloe vera, and Nyctanthes arbor-tristis were compared, as shown in Fig. 3. The red bars indicate the compounds with the highest total IFIE from each plant. It should be noted that anthocyanin (N8) was excluded from the list of Mpro inhibitors because it had a positive charge and interacted differently with the Mpro residues compared with that of the other compounds.
As described in our previous study [23], of the 12 compounds from Moringa oleifera, niaziminin (M9) was found to bind more strongly to the Mpro (Table S4 of SM). This compound interacted strongly with the negatively charged Glu166 residue, resulting in a strong binding to the Mpro. The total IFIE was evaluated, having a magnitude of 136.5 kcal/mol, at least 17.5 kcal/mol higher than that of the other compounds extracted from Moringa oleifera. Therefore, we proposed novel compounds based on niaziminin (M9) to investigate whether the introduction of a hydroxyl group into niaziminin (M9) would enhance its binding affinity to the Mpro. Niaziminin (M9) is a thiocarbamate isolated from the leaves of Moringa oleifera [51]. The presence of an acetoxy group at the 4′-position of niaziminin (M9) is thought to be important for its ability to inhibit the activation of Epstein-Barr virus [52].
As indicated in Fig. 3, of the 34 remaining compounds after the exclusion of N8, compound feralolide (A14) from Aloe vera had the highest total IFIE (150.1 kcal/mol). This IFIE was approximately 13.6 kcal/mol greater than that of Table 2 This table shows the compounds extracted from Nyctanthes arbor-tristis with their lowest binding energy (BE) (kcal/mol), the number of poses, and the Mpro residues involved in H-bonds for the selected clusters obtained using the AutoDock4.2.6 program [40]. The 256 poses created were classified into several clusters based on the structural similarity, and each cluster was ranked in the order of the BE between the Mpro and the compound. We selected three clusters with the highest number of poses and evaluated the total IFIE (kcal/mol) of each compound with all the Mpro residues using the FMO method. These values are listed in the last column To elucidate the reason for the strong binding affinity between feralolide (A14) and the Mpro, we investigated the IFIEs between the Mpro residues and feralolide (A14), eupatorin (A13), and aloesaponarin-I (A12), as these compounds had similar sizes among the 14 compounds extracted from Aloe vera as shown in Fig. 1. Feralolide (A14) (Fig. 4a) bound most strongly to Glu166 and also interacted with other Mpro residues Leu27, Asn142, and Leu167, resulting in a total IFIE of high magnitude. Therefore, these residues were found to be the main contributors to the strong interaction between the Mpro and feralolide (A14). In contrast, aloesaponarin-I (A12) (Fig. 4c) interacted strongly with only Glu166 and weakly with Asn142 and Gln189, resulting in a total IFIE of lower magnitude with the Mpro. As shown in Fig. 4b, the IFIEs of eupatorin (A13) with the Mpro were noticeably different from those of feralolide (A14) and aloesaponarin-I (A12). There was no strong interaction between eupatorin (A13) and the Mpro residues, and only Cys145 interacted with eupatorin (A13) at a higher magnitude than 10 kcal/mol. Therefore, the total IFIE of eupatorine (A13) was of a significantly lower magnitude than that of feralolide (A14) and aloesaponarin-I (A12). The comparison of Figs. 4a-c revealed that Glu166, Leu167, Leu27, and Asn142 were important for the strong binding of feralolide (A14) to the Mpro. Therefore, we investigated the nature of the interactions between the Mpro residues and feralolide (A14), eupatorin (A13), aloesaponarin-I (A12) to understand which functional groups of the compound were responsible for the binding to the Mpro. As shown in Fig. 5a, two hydroxyl groups in the upper end of feralolide (A14) formed strong hydrogen bonds (1.51 and 1.59 Å) with Glu166 and the peptide backbone between Glu166 and Leu167. In addition, the carbonyl group of feralolide (A14) formed a hydrogen bond with the amide group of Asn142, while the hydroxyl group formed a hydrogen bond with the peptide backbone between Thr26 and Leu27. Aloesaponarin-I (A12) (Fig. 5c) formed only one strong hydrogen bond (1.57 Å) with Glu166. In contrast, as shown in Fig. 5b, there was no hydrogen-bonding interaction between eupatorin (A13) and the Mpro residues, although there were only weak electrostatic interactions with Thr26, Gly143, and Cys145. Conversely, Fig. 5a indicates that feralolide (A14) possessed many groups for strong binding with the Mpro residues. In particular, the two hydroxyl groups in the upper end effectively strengthened the interaction between feralolide (A14) and the Mpro residues. Therefore, feralolide (A14) was expected to be a drug lead compound for the development of potent inhibitors against COVID-19.
Aloesaponarin-I (A12) and aloesaponarin-II (A8), shown in Fig. 1, are anthraquinones isolated from the roots of Aloe vera, and their antiviral activity against the human influenza Fig. 3 The total inter-fragment interaction energies (IFIEs) of the compounds with the Mpro residues were evaluated using the ab initio FMO calculations; the results for 12, 14, and 9 compounds from Moringa oleifera [23], Aloe vera, and Nyctanthes arbor-tristis, respectively. The red bars indicate the compounds with the total IFIEs of the highest magnitude among the compounds from each plant; anthocyanin (N8) was excluded from the list of candidate inhibitors because it was positively charged and its interactions with the Mpro residues were purely electrostatic and not comparable with that of the other compounds Fig. 4 The IFIEs between the Mpro residues and the compounds from Aloe vera: (a) feralolide (A14), (b) eupatorin (A13), and (c) aloesaponarin-I (A12). The total IFIEs of these compounds with the Mpro residues are also shown. The red bars indicate the Mpro residues with the negative IFIEs whose magnitude was higher than 10 kcal/mol  virus was demonstrated in a previous study [54]. In addition, Parvez et al. [55] tested the anti-hepatitis B virus (HBV) potential of Aloe vera extract and its anthraquinones in hepatoma cells. They reported that the effect of aloe-emodin (A10) is comparable with that of lamivudine (a nucleoside reverse transcriptase inhibitor approved for the treatment of HBV) and appears to be the most promising natural anti-HBV drug, with CYP3A4 activation enhancing its therapeutic efficacy. Our present results (Fig. 3) indicated that aloesaponarin-I (A12) had the second-highest total IFIE of the compounds from Aloe vera, while aloesaponarin-II (A8) had an average-sized total IFIE. Therefore, it was expected from the total IFIE evaluated by the FMO method that aloesaponarin-I (A12) could also be a potent inhibitor of the Mpro.
In addition, we analyzed the interactions of the other compounds with a size of total IFIE larger than 100 kcal/ mol with the Mpro. As shown in Fig. 3, the Aloe vera compounds caffeic acid (A2), esculetin (A6), aloe-emodin (A10), aloesaponarin-I (A12), and feralolide (A14) came into this category. Since esculetin (A6), a coumarin derivative found in various medicinal plants, inhibits proliferation and induces apoptosis in many types of human cancer cells, it is a promising chemotherapeutic agent [56]. In addition, esculetin (A6) has been widely used as an anti-inflammatory, antioxidant, antibacterial, and anti-diabetic therapeutic [57]. Recently, the anti-HBV activity of esculetin (A6) was investigated, and it was shown to inhibit HBV replication effectively both in vitro and in vivo [58]. Here, of the 14 compounds from Aloe vera, esculetin (A6) exhibited the IFIE of the highest magnitude (73 kcal/mol) with Glu166 and formed two hydrogen bonds. These results indicated that hydrogen bonding to Glu166 is a key for Aloe vera compounds to bind strongly to the Mpro. According to Juárez-Saldívar et al. [59], Glu166 is important for maintaining the shape of the Mpro active site. In addition to this hydrogen bond, most of the Aloe vera compounds participate in electrostatic interactions with Gly143, Ser144, and Cys145. In contrast, Nguyen et al. [60] found that Gly143 of the Mpro is the most important residue for hydrogen bonding with ligands, followed by Glu166, Cys145, and His163. Table 2 lists the results of the docking and FMO calculations for the nine compounds extracted from Nyctanthes arbor-tristis. From the two or three clusters for each compound, we selected the cluster with the total IFIE of the highest magnitude and compared the magnitude of this IFIE for each compound. As shown in Fig. 3, the magnitudes of the total IFIEs of the Nyctanthes arbor-tristis compounds, excluding N8, were significantly lower than those of the Aloe vera compounds. Actually, N8 had the highest total IFIE of all the compounds investigated here. However, this strong interaction came from the electrostatic interaction of the positive charge on N8 with the negatively charged residues of the Fig. 7 The IFIEs of the selected compounds from Nyctanthes arbor-tristis: (a) apigenin (N9) and (b) myristic acid (N3). The total IFIEs of these compounds with the Mpro residues are also shown. The red bars indicate the Mpro residues with the negative IFIEs whose magnitude was higher than 10 kcal/mol Mpro, as shown in Fig. 6. Therefore, the electrostatic nature of this interaction of N8 with the Mpro was significantly different from that of the other compounds; thus, in the present study, we eliminated N8 from the list of Mpro inhibitors.
Of the eight remaining compounds from Nyctanthes arbortristis, apigenin (N9) exhibited a total IFIE of high magnitude (94.8 kcal/mol). As indicated in Fig. 7a, apigenin (N9) interacted strongly only with Glu166 via electrostatic interactions, suggesting that Glu166 is a key residue for stabilizing bound apigenin (N9). However, no residue other than Glu166 interacted strongly with apigenin (N9), resulting in the lower magnitude of its total IFIE compared with that of feralolide (A14) and aloesaponarin-I (A12). Previously, the antiviral activity of the naturally occurring flavonoid apigenin (N9) was investigated in vitro and in vivo against several viruses including enterovirus 71, hepatitis C virus, human immunodeficiency virus, and adenoviruses [61]. Therefore, it was expected that apigenin (N9) could also be an inhibitor of the Mpro, although its binding affinity to the Mpro is likely to be lower than those of feralolide (A14) and aloesaponarin-I (A12). Fig. 9 Visualization of the clustering of compounds using a visualized cluster analysis of protein-ligand interaction [62]. Each compound is represented in the vertical axis. The horizontal axis on the right represents the Mpro residues existing within 10 Å of the compound. Red and blue colors represent attractive and repulsive interac-tions, respectively, and the shading represents the magnitude of the absolute IFIE value. The tree on the left indicates the clustering of compounds based on the IFIEs with the Mpro residues. Compounds in the upper group have strong attractive interaction with Glu166 The IFIEs between Mpro residues and long-chain fatty acids such as pelargonic acid (N2) and myristic acid (N3) were also compared. As indicated in Fig. 7b, myristic acid (N3) interacted strongly only with Glu166 and Leu141 of the Mpro and had a total IFIE of lower magnitude compared with that of apigenin (N9). Therefore, we considered that these long-chain fatty acids are not suitable for the potent inhibitors against the Mpro.
To explain the weaker binding of apigenin (N9) to the Mpro compared with that of feralolide (A14), we analyzed the binding interactions of these compounds with the Mpro residues in the binding site. As shown in Fig. 5a, three hydroxyl groups of feralolide (A14) formed strong hydrogen bonds with the Mpro residues, resulting in the total IFIE of the highest magnitude. In contrast, as indicated in Fig. 8, only one hydroxyl group of apigenin (N9) formed a hydrogen bond with Glu166, although apigenin (N9) has a similar structure to feralolide (A14) and aloesaponarin-I (A12). As a result, the magnitude of the total IFIE of apigenin (N9) was significantly lower than that of feralolide (A14) and aloesaponarin-I (A12), as shown in Fig. 3. These results demonstrated that the number and the position of hydroxyl groups in these compounds significantly affect their specific interactions with the Mpro residues in the binding site.

Specific interactions of the compounds with the main protease
To elucidate the difference in the specific interactions of the 35 natural compounds contained in Moringa oleifera, Aloe vera, and Nyctanthes arbor-tristis with the Mpro residues, we classified these compounds based on their IFIEs using the visualized cluster analysis of protein-ligand interactions [62]. In this analysis, the ligands were classified into different groups depending on the similarity of the IFIEs of each ligand with the Mpro residues in the ligand-binding site. We eliminated N8 as a candidate inhibitor of the Mpro because it is positively charged and its way of interaction with the Mpro residues was significantly different from that of the other compounds.
As shown in Fig. 9, the compounds were classified into two main groups depending on the nature of their interaction with Glu166 in the Mpro binding site. Feralolide (A14), which had the highest total IFIE among the compounds investigated here, was classified in the upper group in Fig. 9, in which all compounds interacted strongly with Glu166. Within this group, only feralolide (A14) formed a subgroup, while the other compounds formed a completely different subgroup, indicating that no other compounds had similar interactions with the Mpro residues compared to feralolide (A14). In contrast, niaziminin (M9), which had the total IFIE of the second-highest magnitude, belonged to the same subgroup of aloesaponarin-I (A12). Therefore, it was revealed from Fig. 9 that the IFIEs between Mpro residues and niaziminin (M9) and aloesaponarin-I (A12) were similar although these compounds had different chemical structures.
To understand the reason for these differences in IFIEs, we compared the interactions of these compounds with the key residues in the Mpro binding site. As shown in Fig. 5a, the hydroxyl groups at one end of feralolide (A14) formed hydrogen bonds with the negatively charged Glu166, whereas that at the other end formed hydrogen bonds with the peptide backbone between Thr26 and Leu27. These hydrogen bonds of feralolide (A14) with two separate sites in the Mpro-binding pocket were achieved by its vertical conformation in the binding pocket, as shown in Fig. 10a. Notably, of all compounds investigated in the present study, only feralolide (A14) was bound in this conformation. Feralolide (A14) also formed a weak hydrogen bond with Asn142. In contrast, Fig. 10b indicates that niaziminin (M9) and aloesaponarin-I (A12) were bound in a perpendicular conformation in the binding pocket, interacting with Glu166, Gln189, and Asn142. As these compounds did not interact with Thr26 and Leu27, their binding properties were not similar to those of feralolide (A14). The magnitude of the total IFIE of aloesaponarin-I (A12) with the Mpro was the third highest among the compounds employed in this study. Therefore, aloesaponarin-I (A12) is also expected to be a potent inhibitor against the Mpro in addition to feralolide (A14) and niaziminin (M9).

Proposal of novel compounds based on feralolide (A14) and niazirinin (M12)
In our previous molecular simulations [23], we investigated the binding properties of 12 compounds extracted from Moringa oleifera with the Mpro and found that, niaziminin (M9) is bound most strongly to the Mpro. The chemical structures of the 12 compounds and the results of the simulations are shown in Fig. S1 and Tables S1 and S2 of SM. In addition, we proposed novel compounds based on niaziminin (M9) and revealed that the introduction of a hydroxyl group into niaziminin (M9) enhanced its binding to the Mpro.
In the present study, we first attempted to propose candidate potent inhibitors of the Mpro based on feralolide (A14), which had the total IFIE of the highest magnitude among the compounds studied here, as shown in Fig. 3. We replaced each hydrogen atom of feralolide (A14) with a hydroxyl group to propose nine novel compounds and investigated their binding properties with the Mpro in the same conformation of feralolide (A14). To obtain the stable structures of the Mpro complexes with these proposed compounds, we optimized their structures using the MM method in water. In addition, the total IFIEs between the Mpro and the compounds were evaluated using the FMO method to propose novel inhibitors with higher binding affinity to the Mpro. As shown in Table 3, the total IFIE was changed significantly depending on the hydroxylation site. In particular, the magnitude of the total IFIEs of feralolide derivatives (A14b), (A14d), and (A14i) were higher than 160 kcal/mol and at least 10 kcal/mol higher than that of feralolide (A14), indicating that the introduction of a hydroxyl group could significantly enhance the binding of feralolide derivatives to the Mpro.
To increase our understanding of this enhancement, we compared the IFIEs of each compound with the Mpro residues. As shown in Fig. 11a, feralolide (A14) interacted strongly with Glu166, Leu167, Leu27, and Asn142. By introducing a hydroxyl group in the d-site of feralolide (A14d), the IFIE between Glu166 and feralolide (A14) was further enhanced by 30 kcal/mol as shown in Fig. 11c, resulting in the total IFIE of A14d having the highest magnitude of 171.1 kcal/mol. In contrast, as shown in Fig. 11b-d, the introduction of a hydroxyl group in the b-site or i-site of feralolide (A14) enhanced the IFIE with Gln189 or Cys145 although the effect was not as strong as that seen with feralolide derivative (A14d). Therefore, the hydroxylation in the d-site of A14 was found to be most effective for enhancing the binding between feralolide (A14) and the Mpro. To understand this enhancement, we investigated the interactions of the proposed feralolide derivatives with some key residues of the Mpro. As shown in Fig. 5a, feralolide (A14) formed hydrogen bonds with Glu166, Asn142, and the peptide backbones between Glu166 and Leu167 and between Thr26 and Leu27. In the complex with feralolide derivative (A14d) (Fig. 12b), the extra hydroxyl group created an additional hydrogen bond with Glu166. As a result, the IFIE of Glu166 was significantly enhanced as shown in Fig. 11c. In contrast, Fig. 12a indicates that if the b-site was hydroxylated to produce feralolide derivative (A14b), this formed a new hydrogen bond with Gln189, leading to the enhancement of IFIE with Gln189. As the hydroxyl group of feralolide derivative (A14i) interacted electrostatically with Cys145, as shown in Fig. 12c, the IFIE of Cys145 was enhanced in feralolide derivative (A14i), as indicated in Fig. 11d. As mentioned above, depending on the site of hydroxylation on feralolide (A14), the interactions of the feralolide (A14) derivatives with some key residues of the Mpro were significantly changed. Of the proposed feralolide (A14) derivatives, the derivative (A14d) was found to have the total IFIE of the highest magnitude with the Mpro residues, indicating that it could be a candidate compound as a potent inhibitor of the Mpro.
Next, we proposed novel derivatives of niazirinin (M12) and investigated their binding properties to the Mpro because niazirinin (M12) had the total IFIE of the secondhighest magnitude of the compounds extracted from Moringa oleifera, as shown in Fig. 3. One hydrogen atom of niazirinin (M12) was replaced by a hydroxyl group to create candidate compounds, and their total IFIEs with the Mpro were investigated. As indicated in Table S5 of SM, only niazirinin derivatives (M12c) and (M12d) had higher total IFIEs compared with niazirinin (M12). However, the magnitude of the total IFIEs of the niazirinin derivatives (M12c) and (M12d) was lower than that (150 kcal/mol) of feralolide (A14). Therefore, it was revealed that the derivatization of niazirinin (M12) by the introduction of a hydroxyl group was not as effective for enhancing the binding affinity to the Mpro, compared with the derivatives of feralolide (A14).
To understand the change in the total IFIE is due to the introduction of an extra hydroxyl group into niazirinin (M12), we compared the IFIEs to the Mpro residues for niazirinin (M12) and its derivatives (M12c), (M12d). As shown in Fig. S2a of SM, niazirinin (M12) interacted strongly with Glu166 and Asn142. The effect of the hydroxylation at the c-site of M12 was, at most, 5 kcal/mol for Cys145 (Fig. S2b). In contrast, the hydroxylation at the d-site significantly enhanced the interactions of niazirinin (M12) with the residues Leu141 and Glu166 (Fig. S2c). However, the interactions of niazirinin (M12) with some residues around Leu141 were weakened by the hydroxylation, resulting only a small enhancement of the total IFIE for the niazirinin derivative (M12d).
Therefore, we investigated the nature of the interactions of Mpro key residues with niazirinin (M12) and derivatives (M12c) and (M12d). As shown in Fig. S3a of SM, niazirinin (M12) formed a strong hydrogen bond with Glu166 and interacted electrostatically with Asn142 and Gly143. As shown in Fig. S3b, the hydroxyl group Fig. 11 (a) The IFIEs of feralolide (A14) with the Mpro residues and the difference in IFIEs between (b) feralolide (A14) and its derivative (A14b), (c) feralolide (A14) and its derivative (A14d), and (d) feralolide (A14) and its derivative (A14i). The red bars indicate the residues with an IFIE difference, the magnitude of which was higher than 10 kcal/mol. These residues interacted more strongly with feralolide derivatives (A14b), (A14d), and (A14i) than feralolide (A14) introduced at the c-site of niazirinin (M12) interacted electrostatically with Cys145. In addition, a hydrogen bond was formed with a distance of 1.65 Å between the Glu166 side chain and a hydroxyl group of another site on the niazirinin derivative (M12c). The electrostatic interactions with the Leu141-Asn142-Gly143 peptide were also enhanced, as shown in Fig. S3b. These changes in interactions in the binding site were considered to cause a significant (about 20 kcal/ mol) increase in the total IFIE of the niazirinin derivative (M12c). In contrast, the hydroxyl group introduced at the d-site of niazirinin (M12d) formed a strong hydrogen bond with the peptide backbone between Phe140 and Leu141 (Fig. S3c). In addition, an extra hydrogen bond was formed between the niazirinin derivative (M12d) and Glu166, resulting in the enhancement of the IFIE for Glu166. However, as shown in Fig. S2b, c, the effect of introducing a hydroxyl group to niazirinin (M12) on the interactions between M12 and the Mpro was found to be at most 12 kcal/mol, which Fig. 12 Models of the molecular interactions of the feralolide derivatives (ball-and-stick models) with the key residues of Mpro (stick models): (a) feralolide derivative (A14b), (b) feralolide derivative (A14d), and (c) feralolide derivative (A14i). The red and blue lines indicate the hydrogen bonding and electrostatic interactions, respectively. The site of hydroxylation in the derivative is marked by a blue or red circle was less significant than that (at most 30 kcal/mol) for the hydroxylated feralolide (A14) derivatives shown in Fig. 11.
To elucidate which compound proposed in our study was bound more strongly to the Mpro, we compared their total IFIEs in Fig. 13. Niaziminin derivative (M9d) and feralolide derivatives (A14d), (A14b) had total IFIEs with a rather higher magnitude compared to those of the other compounds, indicating that they were expected to be the potent inhibitors of the Mpro. To confirm the stability of binding for these compounds to the Mpro, we conducted 100 ns MD simulations for the complexes of Mpro with each of these compounds at 310 K, within an explicit water box using the MD simulation program package GROMACS 2018.4 [63]. The root mean square deviation (RMSD) from the initial structure of the MD simulations were compared in Fig. 14 for the Mpro complexes with A14, A14b, and A14d.
Although the RMSDs increased rapidly to 2-4 Å during the initial 1 ns, there is no significant increase during the MD simulations from 20 to 100 ns. It was confirmed that A14, A14b, and A14d stayed in the ligand-binding pocket of the Mpro even at 100 ns. Figure 14 also indicated that RMSDs for the A14b and A14d complexes were significantly smaller than that of the A14 complex. This result can be explained by the trend in the magnitude of the total IFIEs of the complexes shown in Table 3. Their total IFIEs evaluated in the present study are -171.1 (A14d), -168.7 (A14b), and -150.1 (A14), respectively, indicating that A14d and A14b can bind more strongly to the ligand-binding pocket of the Mpro and stay at the same position compared with A14. In addition, we carried out the same MD simulations for the Mpro complexes with M9/M9d compound. As shown in Fig. S4, the RMSD values of the complexes have no significant increase during the MD simulations from 10 to 100 ns.
Moreover, we investigated the root mean square fluctuation (RMSF) of each Mpro residue and the compound for the Mpro complexes, to elucidate which residues or compounds fluctuate more significantly during the MD simulations. The results for the A14 and M9 derivatives are shown in Figs. 15 and S5, respectively. As indicated in Fig. 15, for the Mpro complexes with the A14 derivatives, although the C-and N-terminal residues had rather large RMSFs, the RMSFs of the other residues are smaller than 4 Å. And the RMSFs of the compounds were found to be 1.46 (A14), 2.61 (A14b), and 0.85 Å (A14d), respectively. Similar results were obtained for the Mpro complexes with the M9/M9d compound as shown in Fig. S5. These MD results revealed that the compounds can stay at the binding pocket of the Mpro even at 310 K. Therefore, we concluded that our proposed A14b, A14d, and M9d compounds are expected to be potent inhibitors against the Mpro. Fig. 13 The total IFIEs of our proposed compounds based on feralolide (A14), niaziminin (M9) [23], and niazirinin (M12) with the Mpro evaluated using ab initio FMO calculations. The red bar indicates the compound with the total IFIE of the highest magnitude among the proposed compounds Fig. 14 Root mean square deviation (RMSD) from the initial structure of the MD simulation for the Mpro complexes with A14, A14b, or A14d compound. The blue, orange, and gray lines represent the RMSD for the complexes with A14, A14b, or A14d, respectively

Conclusions
To find novel potent inhibitors of the SARS-CoV-2 Mpro, we investigated the binding properties between the Mpro and the 35 compounds extracted from Moringa oleifera, Aloe vera, and Nyctanthes arbor-tristis, using molecular simulations based on protein-ligand docking, MM optimizations, and ab initio FMO calculations. The FMO results revealed Fig. 15 Root mean square fluctuation (RMSF) of each Mpro residue and the A14 derivative in the Mpro complexes with (a) A14, (b) A14b, or (c) A14d, respectively. RMSF of the A14 derivative is indicated by the red line at the right end that feralolide (A14) from Aloe vera bound most strongly to the Mpro, making it a promising Mpro inhibitor candidate. To enhance the binding affinity of feralolide (A14), we sequentially substituted a hydroxyl group for each of the hydrogens in the feralolide (A14) structure. We repeated this process for niaziminin (M9) and niazirinin (M12) from Moringa oleifera and investigated the binding properties of their derivatives with the Mpro. The results elucidated that feralolide derivative (A14d) shown in Table 3 as well as our previously proposed niaziminin derivative (M9d) [23] bound strongly to the Mpro. Therefore, these derivatives are expected to be potent inhibitors of the Mpro.