The Inhibitory Potential of Ferulic Acid Derivatives against the SARS-CoV-2 Main Protease: Molecular Docking, Molecular Dynamics, and ADMET Evaluation

The main protease (Mpro) of SARS-CoV-2 is an appealing target for the development of antiviral compounds, due to its critical role in the viral life cycle and its high conservation among different coronaviruses and the continuously emerging mutants of SARS-CoV-2. Ferulic acid (FA) is a phytochemical with several health benefits that is abundant in plant biomass and has been used as a basis for the enzymatic or chemical synthesis of derivatives with improved properties, including antiviral activity against a range of viruses. This study tested 54 reported FA derivatives for their inhibitory potential against Mpro by in silico simulations. Molecular docking was performed using Autodock Vina, resulting in comparable or better binding affinities for 14 compounds compared to the known inhibitors N3 and GC376. ADMET analysis showed limited bioavailability but significantly improved the solubility for the enzymatically synthesized hits while better bioavailability and druglikeness properties but higher toxicity were observed for the chemically synthesized ones. MD simulations confirmed the stability of the complexes of the most promising compounds with Mpro, highlighting FA rutinoside and compound e27 as the best candidates from each derivative category.


Introduction
During the past two years, public health and socioeconomic life have been going through a severe crisis due to the Coronavirus disease 2019 (COVID-19) pandemic caused by severe acute respiratory syndrome Coronavirus-2 (SARS-CoV-2), which has led to more than 450 million cases and more than 6 million deaths worldwide [1]. SARS-CoV-2 belongs to the family of coronaviruses, which have posed a threat to public health in the past, with the severe acute respiratory syndrome-coronavirus (SARS-CoV) and Middle East respiratory syndrome (MERS-CoV) outbreaks in 2002 and 2012, respectively [2].
The development of vaccines was remarkably quick, allowing many countries to initiate the vaccination process in the beginning of 2021 and, therefore, provide a valuable weapon to boost immunity against the virus [3]. As far as drugs and other non-vaccine therapeutic options are concerned, remdesivir, an RNA-dependent RNA polymerase inhibitor, is the only one to receive FDA approval for use in COVID-19 patients [4]. Currently, 15 other products have received emergency use authorizations by the FDA, including the protease inhibitor Paxlovid [5].
However, new variants continue to emerge, affecting the transmissibility of the virus, the impact of the disease, and the immunity against it. Currently, variants beta, gamma, delta, and omicron are labeled as variants of concern (VOCs) by WHO [6]. Among them, the omicron variant is the most capable of antigenic escape, causing concern over the efficacy The role of M pro is the proteolytic cleavage of two viral polyproteins pp1a and pp1ab towards the formation of nonstructural proteins required for further viral reproduction, including its self-release from the aforementioned polyproteins [47,48]. There are more than 11 proteolytic sites, characterized by the sequence (Leu-Gln)-(Ser/Ala/Gly), where the peptide bond being hydrolyzed is the one after Gln [49]. The enzyme functions in a homodimeric form, through a mechanism of nucleophilic addition facilitated by a cysteine-histidine catalytic dyad (Cys145-His41). Each monomer is composed of 306 residues arranged in a polypeptide chain with 3 distinct domains (domain I: residues 8-101; domain II: residues 102-184; and domain III: residues 201-303) ( Figure 1). The active site of the enzyme is a cavity formed between domains I and II, which are made of antiparallel β-barrels. Domain III is composed of five α-helices and contributes to the formation of the dimer [50][51][52][53]. Apart from the presence of catalytic residues His41 and Cys145, it is important to mention the existence of four main subsites forming the active site of M pro , labeled as S1, S1′, S2, and S4. According to Stoddard et al. [2], the S1 subsite consists of the side chains of Phe140, Asn142, Ser144, Cys145, His163, Glu166, His172, and the backbone of Leu141, Gly143, His164, and Met165; the S1′ subsite is formed by the side chains of Thr25, His41, Val42, Asn119, Gly143, Cys145, and the backbone of Thr26; S2 is created by the side chains of His41, Met49, Tyr54, Asp187, and the backbone of Arg 188; and S4 is made up of the side chains of Met165, Leu167, Pro 168, Ala191, Gln192, and the backbones of Glu166, Arg188, and Thr190. Among them, His41, Gly143, Ser144, Cys145, and Glu166 have been pointed out as residues playing a major role in protein-ligand interactions in molecular dynamics studies [14]. Moreover, residues Gly143, Ser144, and Cys145 form an oxyanion hole capable of stabilizing the negative charge of ligands, such as that of the carbonyl oxygen of the scissile bond in the natural substrate [9,10,53].
Regarding the catalytic mechanism of the enzyme, it is suggested that it is established on the formation of a nucleophilic ion pair, through a proton transfer from the thiol group Apart from the presence of catalytic residues His41 and Cys145, it is important to mention the existence of four main subsites forming the active site of M pro , labeled as S1, S1 , S2, and S4. According to Stoddard et al. [2], the S1 subsite consists of the side chains of Phe140, Asn142, Ser144, Cys145, His163, Glu166, His172, and the backbone of Leu141, Gly143, His164, and Met165; the S1 subsite is formed by the side chains of Thr25, His41, Val42, Asn119, Gly143, Cys145, and the backbone of Thr26; S2 is created by the side chains of His41, Met49, Tyr54, Asp187, and the backbone of Arg 188; and S4 is made up of the side chains of Met165, Leu167, Pro 168, Ala191, Gln192, and the backbones of Glu166, Arg188, and Thr190. Among them, His41, Gly143, Ser144, Cys145, and Glu166 have been pointed out as residues playing a major role in protein-ligand interactions in molecular dynamics studies [14]. Moreover, residues Gly143, Ser144, and Cys145 form an oxyanion hole capable of stabilizing the negative charge of ligands, such as that of the carbonyl oxygen of the scissile bond in the natural substrate [9,10,53].
Regarding the catalytic mechanism of the enzyme, it is suggested that it is established on the formation of a nucleophilic ion pair, through a proton transfer from the thiol group of Cys145 to the imidazole of His41. The catalytic cysteine attacks the carbonyl of the scissile bond, leading to a thiohemiketal intermediate, while the protonated histidine attacks the N-atom of the peptide bond, creating the acyl-enzyme complex intermediate. The participation of a water molecule in the reaction is of great importance, as it attacks the carbonyl carbon of the substrate's Gln while the catalytic His is being reprotonated, and also stabilizes the polar contacts between residues His41, His164, and Asp187 by interacting with them. The last step of the mechanism is described by the release of Cys145 through the cleavage of its covalent bond with the peptide [53].

Docking Validation
Both known inhibitors examined, N3 and GC376, led to a docking output comparable to the co-crystallization structure of their complex with M pro . N3 is the most widely accepted and analyzed inhibitor of M pro in the literature [51,54,55]. It is often used as a positive control to provide some reference values with which the binding energy and interactions of a purported inhibitor with M pro can be compared. The binding energy for N3 calculated in this work by Vina is −8.26 kcal/mol, whereas Das et al. [56] reported it to be −7.7 kcal/mol and Ahmed et al. [57] −7.5 kcal/mol, a difference that can be attributed to the different pretreatment of the receptor and the ligand structures before the simulation. Superimposing of the co-crystallized protein-ligand complex and the one calculated by the simulation reveals high similarity between the conformation of N3 in the binding site in the two cases, with an RMSD of 3.16 Å (Figure 2a). In the docking output, the interaction of N3 with Cys145 is a 2 Å hydrogen bond, formed between the pentacyclic ring of N3 and the hydrogen attached to the sulfur atom of the protein residue. Whereas, as determined by the co-crystallization data, it is a 1.8 Å covalent bond between the sulfur atoms of Cys145 and the Cβ atom of the vinyl group of the inhibitor. Although a different part of the inhibitor is bound to the protease, it results in a very similar geometry of the molecule. Moreover, among the seven protein-ligand hydrogen bonds calculated for the co-crystallization structure, involving residues Gly143, His163, His164, Glu166, Gln189, and Thr190, three were also calculated from the simulation (including Glu166 and Gln189). GC376 is another broad-spectrum inhibitor that has demonstrated activity against the main protease of various coronaviruses [58]. In the present simulation, the binding energy calculated for GC376 was −7.80 kcal/mol. The docking output is remarkably similar to the co-crystallized structure, with the ligands in the two conformations having an RMSD as low as 1.0 Å (Figure 2b). In addition, another study featuring the docking of GC376 to PDB 6LU7 with AutoDock Vina reported a binding energy of −8.1 kcal/mol [59], which is very close to the result of our study. The hydrogen bonds that formed appear to be matching to a great extent, as GC376 interacts with His41, Phe140, His164, Glu166, and Gln189 in the co-crystallized complex and with His41, Phe140, His163, Glu166, and Gln189 in the docked complex.
Biomedicines 2022, 10, x FOR PEER REVIEW 7 of 45 Figure 2. Superimposed binding modes of the inhibitors N3 (a) and GC376 (b) as occurring from the co-crystallization structure (PDB ID: 6LU7 and 7D1M, respectively) (petrol blue) and the molecular docking simulation (purple), respectively. The root-mean-square-deviation (RMSD) of the binding between the co-crystallized and docked complex is indicated.
Overall, the docking results for the known inhibitors are highly comparable to the binding modes observed in reality, reinforcing the validity of molecular docking as a tool to provide indications for the inhibitory potential of the compounds under evaluation.

FA as a 'Reference' Ligand
Although FA derivatives have not been extensively studied in silico for their inhibitory potential against SARS-CoV-2 M pro , FA itself has been part of related studies featuring . Superimposed binding modes of the inhibitors N3 (a) and GC376 (b) as occurring from the co-crystallization structure (PDB ID: 6LU7 and 7D1M, respectively) (petrol blue) and the molecular docking simulation (purple), respectively. The root-mean-square-deviation (RMSD) of the binding between the co-crystallized and docked complex is indicated.
Overall, the docking results for the known inhibitors are highly comparable to the binding modes observed in reality, reinforcing the validity of molecular docking as a tool to provide indications for the inhibitory potential of the compounds under evaluation.

FA as a 'Reference' Ligand
Although FA derivatives have not been extensively studied in silico for their inhibitory potential against SARS-CoV-2 M pro , FA itself has been part of related studies featuring natural compounds. The binding energy calculated by Autodock Vina for FA in this work was −5.88 kcal/mol for the first among six emerging clusters. The next two clusters did not deviate much from this value, with binding energies of −5.60 and −5.47 kcal/mol, respectively ( Figure 3). In the first and third cluster, the phenolic group is situated in the S2 subsite, but the tail of the molecule extends towards the S4 subsite in the first case and towards the S1 subsite in the second. In the case of the second cluster, the phenolic ring is located between the S1 and S1 subsites, in front of the catalytic dyad, while the other end of the molecule is inserted into the S2 subsite. His41 stands out as an interacting residue in all cases, involved in both hydrophobic and pi-pi interactions, while hydrogen bond interactions involve residues closer to the catalytic dyad in the case of clusters 2 and 3, and S4 site residues for cluster 1. Using the protease PDB structure 6LU7 and the Molecular Operating Environment 2019.0102 (MOE) software for molecular docking, two different studies reported almost identical binding energies of −5.33 [60] and −5.35 kcal/mol [61], and the formation of two and one hydrogen bonds, respectively, with residues Thr190 and Glu199, and Glu166. The latter also provided the orientation of the molecule when bound to the active site, which resembles the first cluster of this work (Figure 3), without being identical, however.  These results are quite similar with the results produced in the present work for the first cluster of FA, as the binding energies do not differ significantly and there is one common residue involved in hydrogen bonding (Thr190). Another work reported a binding energy of −10.63 kcal/mol and hydrogen bonds with residues Leu141, Gly143, Ser144, Cys145, His163, and Glu166 [62] while the orientation of FA in the active is identical to the one observed in the second cluster for FA in this work. In accordance with the present results, a docking score of −6.0 kcal/mol was reported for FA when the Autodock Vina and PDB structure 6W63 were used [63]. Overall, the docking scores for FA are quite low compared to established inhibitors such as the previously mentioned N3 (−8.26 kcal/mol) and GC376 (−7.80 kcal/mol). Therefore, the screening of FA derivatives aims to explore their potentially improved properties both in M pro inhibition and in terms of pharmacokinetics.

Docking of Enzymatically Synthesized FA Derivatives
Enzymatic synthesis of FA derivatives has been prevalently dominated by esterification of FA or transesterification of a respective activated donor (such as methyl ferulate, These results are quite similar with the results produced in the present work for the first cluster of FA, as the binding energies do not differ significantly and there is one common residue involved in hydrogen bonding (Thr190). Another work reported a binding energy of −10.63 kcal/mol and hydrogen bonds with residues Leu141, Gly143, Ser144, Cys145, His163, and Glu166 [62] while the orientation of FA in the active is identical to the one observed in the second cluster for FA in this work. In accordance with the present results, a docking score of −6.0 kcal/mol was reported for FA when the Autodock Vina and PDB structure 6W63 were used [63]. Overall, the docking scores for FA are quite low compared to established inhibitors such as the previously mentioned N3 (−8.26 kcal/mol) and GC376 (−7.80 kcal/mol). Therefore, the screening of FA derivatives aims to explore their potentially improved properties both in M pro inhibition and in terms of pharmacokinetics.

Docking of Enzymatically Synthesized FA Derivatives
Enzymatic synthesis of FA derivatives has been prevalently dominated by esterification of FA or transesterification of a respective activated donor (such as methyl ferulate, MFA or vinyl ferulate, VFA), under a low water content towards the production of stabilized esters with tailored lipophilicity [36,64,65]. Esterases that are widely used for this purpose are the triaglycerol lipases (EC 3.1.1.3) due to their broad specificity towards glycerides and related substrates. Nevertheless, feruloyl esterases (EC 3.1.1.73), which have specificity towards hydroxynnamic acids, have also been employed for this purpose as they can offer the unique advantage of catalyzing the synthesis of esters with a wide range of substitutions on the phenolic ring, resulting in either more lipophilic (e.g., alkyl esters) or hydrophilic esters (e.g., sugar esters) [66][67][68][69][70]. Transglycosylation of FA with the glucoside rutin has been described [35]. An overview of the major FA derivatives that have been synthesized enzymatically along with the representative synthesis routes is described in Table 1. These derivatives were selected for further in silico simulations. Table 2 presents the results of the docking simulations for the major enzymatically synthesized derivatives.

Alkyl and Alkenyl Esters of FA
The FA derivatives belonging to this category that were investigated include methyl, ethyl, propyl, butyl, isobutyl, pentyl, isopentyl, prenyl, hexyl, octyl, dodecyl, and octadecyl ferulates. Overall, the calculated binding energies ranged from −5.20 to −6.75 kcal/mol, with prenyl ferulate exhibiting the highest binding affinity to the active site of M pro . Although there is not a definite correlation between the binding energy and the potential of in vitro inhibitory effect, it is not encouraging that these values are significantly lower than the inhibitors N3 and GC376. The total results in terms of the binding energies and interactions calculated are shown in Table 2. It is interesting that prenyl ferulate is the only alkenyl ester among this category, and has better binding compared to the respective alkyl ester (isopentyl ferulate, −6.58 kcal/mol). It is observed that the binding energies become progressively lower as the carbon chain of the substitutions becomes larger, until the point of five carbons. For the larger substitutions, the binding affinity starts to decrease again.
The prevalent conformation observed is the one where the phenolic ring of the FA moiety is stabilized at the S1 subsite, in front of the catalytic residue Cys145, while the carbon chain extends towards the S2/S4 subsites. This particular orientation is seen in the first cluster for all ligands in this category, except for pentyl ferulate, where the orientation corresponds to the second cluster ( Figure 4). The binding energy is almost identical between the two clusters. Moreover, no such conformation is observed in the cases of dodecyl and octadecyl ferulate, which possess a long carbon chain that folds differently into the active site.
Another commonly occurring geometry is characterized by the phenolic moiety of the ligands being located at the S2 subsite while the rest of the molecule is orientated towards the S1 or S1 subsites, as seen, for example, in the case of the third cluster for methyl ferulate or the second cluster for butyl ferulate, respectively ( Figure S1). However, taking into consideration the facts that this conformation does not appear as consistently as the previously described one and, when it does, it represents the second or third clusters emerging from the simulation, which often have a considerably lower binding affinity to the active site of the protease, it can be suggested that the most likely position for the phenolic ring is the one in front of the catalytic cysteine. In this spatial arrangement of the molecule, the hydroxyl substitution, which readily contributes to the formation of hydrogen bonds, appears to play an important role in the stabilization of the molecule in a position where access to one of the two catalytic residues is restricted, therefore potentially resulting in more effective inhibition of the catalytic activity of the enzyme.     carbon chain extends towards the S2/S4 subsites. This particular orientation is seen in the first cluster for all ligands in this category, except for pentyl ferulate, where the orientation corresponds to the second cluster ( Figure 4). The binding energy is almost identical between the two clusters. Moreover, no such conformation is observed in the cases of dodecyl and octadecyl ferulate, which possess a long carbon chain that folds differently into the active site.  Regarding comparable results in the literature, the binding of methyl ferulate to M pro (PDB ID:6W63) was simulated with AutoDock 4.2, which calculated a binding energy of −5.19 kcal/mol, and the formation of hydrogen bonds with the residues Thr190 and Gln192 [100]. The binding energy is comparable to the one calculated in this work (−5.73 kcal/mol), but the hydrogen bond interactions are located in different parts of the active site, as the present study indicated interactions around the catalytic dyad with residues Gly143, Ser144, and Cys145 for the first cluster.

Fatty Esters of FA and Other Related Esters
This category of compounds includes oleyl, glyceryl, diglyceryl, tocopheryl, and βsitosteryl ferulates. The binding energies calculated via the docking simulation range from −5.17 to −7.81 kcal/mol. The lowest binding affinity was predicted for oleyl ferulate, which has an 18-carbon aliphatic chain attached to the FA moiety and exhibits similar results to these of octadecyl ferulate, which possesses an equally long carbon chain, both in terms of the binding energy (−5.17 and −5.20 kcal/mol, respectively) and the conformation of the first cluster. Glyceryl, diglyceryl, and tocopheryl ferulates do not fall below −7.0 kcal/mol, thus exhibiting a lower binding affinity compared to the inhibitors N3 (−8.26 kcal/mol) and GC376 (−7.80 kcal/mol), which is, however, not high enough to exclude the possibility of effective inhibition. The most promising results were yielded for β-sitosteryl ferulate in the first cluster, where the triterpene moiety of the molecule occupies the S4 and S1 subsites while the FA moiety extends to the S1 subsite and out of the active site cavity. This ligand emerges as particularly promising, as its binding energy (−7.81 kcal/mol) is very close to that of GC376 (−7.80 kcal/mol). β-Sitosteryl ferulate also forms three hydrogen bonds, with residues Thr25, Gly143, and the catalytic His41. Another molecular docking study also using Autodock Vina and PDB structure 6LU7 reported an identical binding energy of −7.8 kcal/mol for β-sitosteryl ferulate; however, in this case, the compound appears to bind to M pro at a different site [101].
Regarding other patterns observed in the binding mode to the active site, it is observed that the smaller ligands in this category (glyceryl and diglyceryl ferulate) bind to the active in a similar manner as the FA derivatives described in the previous paragraph, where the FA phenolic moiety is situated in front of Cys145 at the S1 subsite, and its hydroxyl group participates in hydrogen bonding with the neighboring residues while the tail of the molecule is orientated towards the S4/S2 subsites. This orientation is seen in the first cluster of both glyceryl and diglyceryl ferulate. As far as the molecules with the larger carbon chains are concerned (oley and tocopheryl ferulates), no binding patterns were observed while the orientation of tocopheryl ferulate in all three clusters appeared to be more efficiently blocking access to the catalytic dyad compared to oleyl ferulate, where the FA moiety mostly occupied the active site. The structures of the ligands and the conformations described above are shown in Figure 5 for cluster 1 and Figure S2 for clusters 2 and 3.
The range of binding energies calculated for the monosaccharide-based esters is quite narrow, ranging from −7.06 kcal/mol for D-fructose ferulate, which is the only FA derivative with a keto-hexose tested, to −7.37 kcal/mol for D-xylose ferulate. In comparison to the inhibitors N3 and GC376, these results are not exceptional, but they are comparable. Therefore, also taking into consideration the improved solubility of these compounds, these sugar esters of FA are worthy of further investigation. As far as the first clusters are concerned, two orientations are more often seen. One is characterized by the positioning of the penta-or hexacyclic ring in the S1 subsite, close to the catalytic Cys145, and the extension of the rest of the molecule horizontally, towards the S4 (or sometimes S2) subsite, as in the case of the first cluster for D-galactose ferulate ( Figure 6). With the exception of D-glucose and D-mannose ferulates, where the FA moiety occupies the S2 subsite, and L-arabinose ferulate, for which this conformation is not observed at all in the first three clusters, the interactions observed are almost identical and include hydrogen bonds with residues Leu141, Asn142, Gly143, Ser144, Cys145, His163, Arg188, Thr190, and Gln192. The second most prevalent orientation resembles the dominant orientation described in the previous paragraphs, and is opposite to the aforementioned one, with the phenolic ring of FA being stabilized in front of Cys145 at the S1 subsite and the sugar substitution blocking the S4 subsite. The residues more often involved in hydrogen bonding in this case are quite similar, including Leu141, Gly143, Ser144, Glu166, Thr190, and Gln192. Overall, based on the frequency of occurrence of the presented conformations in the first three clusters resulting from the simulation (Figures 6 and S3), the first binding mode described is the prevailing one. This supposition is supported by the fact that in the cases where this geometry appears in the second cluster instead of the first, the resulting binding energies of the two clusters are almost identical. An example is D-glucose ferulate, with a binding energy of −7.09 kcal/mol for the first cluster and −7.08 kcal/mol for the second cluster.
The sugar esters of FA also include coupling with disaccharides D-lactose, D-sucrose, D-maltose, D-cellobiose, and xylobiose. Their binding energies provide a positive indication of the inhibitory potential, as they fluctuate within a small range comparable to that of the reference inhibitors, from −7.47 kcal/mol for D-maltose ferulate to −7.97 kcal/mol for xylobiose ferulate. Three binding patterns stand out from the evaluation of the first three clusters for each compound, with only one of them being present in all the cases. This involves the FA moiety blocking the S4 subsite, the monosaccharide closer to FA taking up the S1 subsite, and the second monosaccharide extending upwards at the S1 subsite, as seen in the case of the third cluster for D-lactose ferulate or the first cluster for D-sucrose ferulate (Figure 7).
The occupation of the same regions of the active site is, however, not translated into identical hydrogen bond interactions. The most commonly occurring ones involve residues Gly143 and Thr190 in four out of five ligands and Thr26 in three out of five. Another conformation observed in four out of the five compounds (e.g., the first cluster for Dmaltose ferulate) is a vertical one, with the cyclic ring of the sugar closer to FA blocking the catalytic cysteine between the S1 and S1 subsites and the FA extending towards S1 . An evident motif is observed in the interactions between the compounds in this orientation and the protease active site, with Ser46, Gly143, Ser 144, Cys145, and Glu166 emerging as hydrogen bond hotspots. Lastly, a conformation involving the occupation of all four subsites, S4 by the phenolic ring of FA, S2 by its carbonyl group, S1 by the monosaccharide ring closer to FA, and S1 by the second monosaccharide, respectively, is observed in two of the five cases (in the second cluster of D-lactose ferulate and the third clusters of D-sucrose ferulate ( Figure S4). This binding pattern is stabilized by hydrogen bonds of the ligands with residues Phe140, Leu141, Asn142, Gly143, His163, and Glu166. substitution blocking the S4 subsite. The residues more often involved in hydrogen bonding in this case are quite similar, including Leu141, Gly143, Ser144, Glu166, Thr190, and Gln192. Overall, based on the frequency of occurrence of the presented conformations in the first three clusters resulting from the simulation (Figure 6, Figure S3), the first binding mode described is the prevailing one. This supposition is supported by the fact that in the cases where this geometry appears in the second cluster instead of the first, the resulting binding energies of the two clusters are almost identical. An example is Dglucose ferulate, with a binding energy of −7.09 kcal/mol for the first cluster and −7.08 kcal/mol for the second cluster. The sugar esters of FA also include coupling with disaccharides D-lactose, D-sucrose, D-maltose, D-cellobiose, and xylobiose. Their binding energies provide a positive indication of the inhibitory potential, as they fluctuate within a small range comparable to that of the reference inhibitors, from −7.47 kcal/mol for D-maltose ferulate to −7.97 kcal/mol for xylobiose ferulate. Three binding patterns stand out from the evaluation of the first three clusters for each compound, with only one of them being present in all the cases. This involves the FA moiety blocking the S4 subsite, the monosaccharide closer to FA taking up the S1 subsite, and the second monosaccharide extending upwards at the S1′ subsite, as seen in the case of the third cluster for D-lactose ferulate or the first cluster for D-sucrose ferulate (Figure 7). Two other disaccharide derivatives include two FA groups, galactobiose ferulate and arabinobiose ferulate, with very promising binding energies of −8.36 and −7.88 kcal/mol, respectively, the former of which is the second highest among the enzymatic derivatives tested. When observing the first three clusters for the two ligands, only the first cluster for galactobiose ferulate had an orientation comparable to any of the binding patterns detected for the other disaccharide derivatives. This conformation involves the two monosaccharides occupying the S1 and S1 subsites and one of the FA groups extending to the S4 subsite while the other is located above the S1 subsite, emerging out of the active site cavity (Figure 7).
The only trisaccharide esterified with one FA moiety tested is raffinose ferulate, which exhibited the third best binding affinity (−8.34 kcal/mol) to M pro among the enzymatic derivatives studied, which also surpasses that of the native inhibitor N3. These results were only exceeded by galactobiose ferulate and one of the structures examined for FOS ferulate (FOS ferulate 3), which will be further described below. Raffinose ferulate forms eight hydrogen bonds with residues Thr26, Ser144, Cys145, His163, Glu66, Thr190, and Gln192, most of which are between the hydroxyl groups of the sugar moieties and residues located at the S1 and S4 subsites. The binding mode of the compound as predicted by the docking simulation appears to be limiting access to Cys145. (Figure 8). Another category of compounds tested is sugar alcohol esters, and, more specifically, the FA esters with D-mannitol, D-sorbitol, and D-xylitol. These compounds had the lowest binding affinity to the active site among the sugar derivatives of FA: −6.38, −6.43, and −6.68 kcal/mol, respectively. In this category, the conformation where the FA ring is in front of Cys145, between the S1 and S1 subsites, while the polyol is orientated towards the S4 subsite is the one that prevails in the first clusters and is stabilized in all three cases with hydrogen bonds with residues Leu141, Gly143, Ser144, Cys145, Glu166, Arg188, Thr190, and Gln192 ( Figure 8). The occupation of the same regions of the active site is, however, not translated into identical hydrogen bond interactions. The most commonly occurring ones involve residues Gly143 and Thr190 in four out of five ligands and Thr26 in three out of five. Another conformation observed in four out of the five compounds (e.g., the first cluster binding affinity to the active site among the sugar derivatives of FA: −6.38, −6.43, and −6.68 kcal/mol, respectively. In this category, the conformation where the FA ring is in front of Cys145, between the S1 and S1′ subsites, while the polyol is orientated towards the S4 subsite is the one that prevails in the first clusters and is stabilized in all three cases with hydrogen bonds with residues Leu141, Gly143, Ser144, Cys145, Glu166, Arg188, Thr190 and Gln192 (Figure 8). Fructooligosaccharides (FOSs) are oligomers of fructose units linked by beta (2-1) glycosidic bonds, with the terminal unit often being a glucose. These oligosaccharides can comprise 2 to 60 monomers [102]; therefore, their exact structure is unknown. In this work Fructooligosaccharides (FOSs) are oligomers of fructose units linked by beta (2-1) glycosidic bonds, with the terminal unit often being a glucose. These oligosaccharides can comprise 2 to 60 monomers [102]; therefore, their exact structure is unknown. In this work, the structures of feruloylated FOS fragments proposed by Couto et al. [99] were adopted for the purposes of the simulation. It is stated in the original work that the position of the FA moieties on the oligosaccharide backbone is not confirmed, and the structures proposed are indicative. They include three forms of FOS ferulate, here referred to as FOS ferulate 1, 2, and 3. FOS ferulate 1 is a fragment that is composed of two FA moieties esterified to the O4 and O5 positions of a fructose and exhibits a binding energy of 8.175 kcal/mol, which is significantly high if compared with the rest of the monosaccharide FA esters. This value is very encouraging, being also highly comparable to the respective values for the reference inhibitors. Although FOS is not a monosaccharide, the orientation of the studied monosaccharide fragment, FOS ferulate 1, in its first cluster resembles the prevalent conformation described above for the rest of the monosaccharide-based derivatives, where the monosaccharide is located below Cys145 and the one FA group occupies the S4 subsite while the other extends upwards to block the S1 subsite ( Figure 9). Hydrogen bonds with residues Leu141, Gly143, Ser144, His163, and Thr190 are also common interactions in both cases. FOS ferulate 2 is a fragment depicting two fructose moieties linked with a (1-2) glycosidic bond, where FA is esterified to the O5 position of the first monomer. It exhibited an encouraging binding energy of −7.34 kcal/mol and 8 hydrogen bonds with residues commonly involved in ligand binding, namely Thr26, Gly143, His164, Glu166, and Thr190. As in the case of FOS ferulate 1, although FOS ferulate 2 is an arbitrarily defined fragment rather than a disaccharide, the binding mode of its first cluster can be compared to that of the disaccharides described above, as it follows the predominant pattern observed, where FA blocks the S4 subsite, the fructose linked to it occupies the S1 subsite, and the other fructose monomer extends to the S1 subsite. FOS ferulate 3 is the largest of the molecules tested, including a fructose trimer with four FA groups esterified to the O5 position of the first monomer, the O6 of the second, and the O4 and O6 of the third, respectively. It showed the highest binding affinity (−8.52 kcal/mol), utilizing its volume to take up the entire active site of M pro and being stabilized through 14 hydrogen bonds, with residues Ser46, Leu141, Asn142, Gly143, Ser144, Cys145, His163, His164, and Gln192.
follows the predominant pattern observed, where FA blocks the S4 subsite, the fructose linked to it occupies the S1 subsite, and the other fructose monomer extends to the S1 subsite. FOS ferulate 3 is the largest of the molecules tested, including a fructose trimer with four FA groups esterified to the O5 position of the first monomer, the O6 of the second, and the O4 and O6 of the third, respectively. It showed the highest binding affinity (−8.52 kcal/mol), utilizing its volume to take up the entire active site of M pro and being stabilized through 14 hydrogen bonds, with residues Ser46, Leu141, Asn142, Gly143 Ser144, Cys145, His163, His164, and Gln192.  Conformations for the second and third cluster results for the simulations regarding trisaccharide, polyol, and polysaccharide esters are presented in Figures S5 and S6.

FA Glycosides
FA rutinoside is a rutinase-synthesized glycoside, which exhibited an exceptionally high binding affinity and a binding energy of 8.404 kcal/mol. It forms several hydrogen bonds with residues Thr24, Thr25, Thr45, Leu141, Gly143, Ser144, Cys145, His163, Glu166, Arg188, and Gln189. The docking simulation output suggests that it binds to the active site of M pro blocking the S1, S1 , and S4 subsites, with its glucosyl group being stabilized through multiple hydrogen bonds at the S1 subsite, the rhamnosyl group occupying the S4 subsite, and the FA moiety extending upwards into the S1 subsite ( Figure 10). The conformations of the second and third cluster are shown in Figure S7. Conformations for the second and third cluster results for the simulations regarding trisaccharide, polyol, and polysaccharide esters are presented in Figures S5 and S6.

FA Glycosides
FA rutinoside is a rutinase-synthesized glycoside, which exhibited an exceptionally high binding affinity and a binding energy of 8.404 kcal/mol. It forms several hydrogen bonds with residues Thr24, Thr25, Thr45, Leu141, Gly143, Ser144, Cys145, His163, Glu166, Arg188, and Gln189. The docking simulation output suggests that it binds to the active site of M pro blocking the S1, S1', and S4 subsites, with its glucosyl group being stabilized through multiple hydrogen bonds at the S1 subsite, the rhamnosyl group occupying the S4 subsite, and the FA moiety extending upwards into the S1′ subsite ( Figure 10). The conformations of the second and third cluster are shown in Figure S7.

Docking of Chemically Synthesized FA Derivatives
The chemically synthesized FA derivatives included in this study were selected based on their promising results in published works regarding antiviral effects ( Table 3). The studies from which the derivatives were selected involve the synthesis of several

Docking of Chemically Synthesized FA Derivatives
The chemically synthesized FA derivatives included in this study were selected based on their promising results in published works regarding antiviral effects ( Table 3). The studies from which the derivatives were selected involve the synthesis of several derivatives and their in vitro evaluation against a range of viruses. The ones selected in this work for investigation are a selection of the compounds that showed the best antiviral potential against the respective tested virus. Such derivatives fall into the categories of FA amides, FA oligomers, and FA derivatives with dithioacetal, acylhydrazone, quinazoline, and chalcone moieties.
The considerable variety in the structure of the ligands in this category is reflected in a broad range of binding energies (from −6.27 to −8.40 kcal/mol) and a wide diversity of binding modes. The binding energy and detailed interactions for each compound are presented in Table 4. A general observation about the binding of the chemically synthesized derivatives, as seen in Figure 11, is that in a considerable number of them, for cluster 1, the phenolic moiety in FA appears to be located between the S2 and S1 subsites, more towards S2. In this position, it is in the close vicinity of the catalytic dyad, as seen, for example, in the cases of compounds 7a, 2, and g18. Particularly in the cases of compound D4, FA dimer, and trimers, the phenolic ring is stabilized closer to the S1 subsite, blocking access to the area around the catalytic residues even more effectively. The binding modes of the second and third clusters were also taken into consideration, but no evident binding pattern was observed ( Figure S8).
The compounds that stood out are two FA derivatives with a quinazoline moiety, namely e27 and e28, with binding energies of −8.11 and −7.76 kcal/mol respectively; 4n, an FA amide and myricetin derivative with a binding energy of −7.82 kcal/mol; F3, an FA-chalcone ester with a binding energy of −7.80 kcal/mol; and triferulic acid, which exhibited a binding energy of −8.32 kcal/mol. The two quinazoline derivatives form hydrogen bonds with His41 and Gln192; the myricetin derivative with residues Leu141 and Gly143; the chalcone derivative with Thr24, Thr 25, Thr45, and Ser46; and the FA trimer with Thr26, Tyr54, and Asp187. Moreover, Gln189 and Met165 often participate in hydrophobic interactions while His41 appears to be involved in pi-pi stacking. Although the interactions do not exhibit any evident pattern, the ligands in all the cases appear to be occupying the area in front of the catalytic dyad, providing a positive indication of potential inhibitory activity. Compounds e27 and e28 have a very similar structure, which results in an almost identical binding mode. It is also interesting that in the case of compound 4n, where the respective substitution is linked to FA at the same position of FA as in the case of e27, and e28 (at C4 of the phenolic ring) binds to the active site in a way that the FA moiety has the same position as in the quinazoline derivatives. Its bulkier substitution, however, allows for higher coverage of the active site and occupation of the S1 subsite as well, as opposed to only S4, S2, and S1 in the case of e27 and e28. The FA trimer also appears to be quite flexible and bends in a way that allows blocking of all four subsites while F3 mostly occupies the S4 and S1 subsites and is stabilized through hydrogen bonds in the upper part of the latter.  FA derivatives with a quinazoline moiety (A) e27
64.80% 28 [105]        The compounds that stood out are two FA derivatives with a quinazoline moiety, namely e27 and e28, with binding energies of −8.11 and −7.76 kcal/mol respectively; 4n, an FA amide and myricetin derivative with a binding energy of −7.82 kcal/mol; F3, an FAchalcone ester with a binding energy of −7.80 kcal/mol; and triferulic acid, which exhibited a binding energy of −8.32 kcal/mol. The two quinazoline derivatives form hydrogen bonds with His41 and Gln192; the myricetin derivative with residues Leu141 and Gly143; the chalcone derivative with Thr24, Thr 25, Thr45, and Ser46; and the FA trimer with Thr26, Tyr54, and Asp187. Moreover, Gln189 and Met165 often participate in hydrophobic interactions while His41 appears to be involved in pi-pi stacking. Although the interactions do not exhibit any evident pattern, the ligands in all the cases appear to be occupying the area in front of the catalytic dyad, providing a positive indication of

Bioavailability and Druglikeness
The estimation of the druglikeness of the compounds under investigation, performed using the SwissADME and MOLSOFT online servers, is presented in Table 5. The size, flexibility, polarity, and solubility/lipophilicity of a compound, and its ability to form hydrogen bonds, affect how easily it can permeate membranes, be absorbed, distributed, and reach and bind to biological targets. More specifically, low lipophilicity hinders the transfer of a compound through cell membranes, but high lipophilicity is usually related to low absorption and high metabolic turnover, and potential toxicity [113]. The majority of the sugar derivatives are less lipophilic than FA while the aliphatic and chemically synthesized derivatives are more lipophilic.
Water solubility is principally a desired property. Almost all the FA derivatives have acceptable solubility, but the derivatization itself causes an increase in the solubility compared to FA mostly in the case of the sugar derivatives. As expected, the alkyl and alkenyl derivatives become progressively less soluble in water as the carbon chain size of the substitution increases. Oleic acid, tocopherol, and sitosterol esters are poorly soluble in aqueous media. As far as the chemically synthesized derivatives are concerned, they display moderate to low solubility. Hydrogen bond donors and acceptors can contribute to the solubility and binding to the desired targets, but the presence of too many of such groups negatively impacts the permeability through cell membranes [114]. The derivatives tested possess an adequate amount of hydrogen bond donors and acceptors, aside from FA sugar derivatives, in which a high number of such groups is observed.
According to Daina et al. [115], the optimum lipophilicity, size, polarity, insolubility, insaturation, and flexibility for a compound to be considered as having good bioavailability correspond to logP o/w values ranging between −0.7 and +5.0, molecular weight between 150 and 500 g/mol, topological polar surface area (TPSA) between 20 and 130 Å 2 , logS between 0 and 6, fraction of carbons in the sp 3 hybridization between 0.25 and 1, and no more than 9 rotatable bonds. Based on these rules, more than half of the compounds have favorable properties in most of the categories. However, the most positive oral bioavailability indications are for the compounds that have less good binding scores.
Another rule for easily evaluating oral bioavailability is the Lipinski's rule of 5, according to which a compound is likely to be orally bioavailable if it has less than 5 hydrogen bond donors and less than 10 acceptors, molecular weight lower than 500 g/mol, and an octanol-water partition coefficient lower than 4.15. It is hopeful that 41 out of the 57 compounds satisfy the Lipinski s rule of 5 regarding oral bioavailability, even though the best hits from the molecular docking simulation are not included in those. However, it is important to note that this is not a strict criterion that can certainly exclude a compound from further investigation. For example, 13.6% of the oral and 50.0% of the non-oral drugs with low solubility and permeability violate the rule [116]. The bioavailability score is another index that demonstrates a similar tendency for the compounds tested as the two previously described criteria.   Overall, the druglikeness scores calculated for the investigated compounds, based on the respective index calculated by MOLSOFT, do not exclude them from being used as drugs. The index is calculated by comparing the molecular properties of the respective compound with a library of selected drug and non-drug compounds [117] and provides an indication of how likely a compound is to function as an oral drug. Scores between 0 and 1 are stronger indicators of druglikeness, but also negative values, particularly over −1, correspond to resemblance to drug molecules. The alkyl and alkenyl derivatives had the most negative scores, the lowest being −0.76, while the FA sugar esters and chemically synthesized derivatives showed more encouraging results, with tocopheryl and sitosteryl ferulate and compound 5 yielding the best scores, with druglikeness scores of 1.14, 1.21, and 1.40, respectively. It is worth mentioning that isobutyl and isopentyl ferulates stood out among the alkyl derivatives in terms of druglikeness, the sugar polyol esters exhibited lower scores compared to the other sugar derivatives, and the chemically designed molecules had an overall higher druglikeness, something that could be expected given the fact that many of them emerged from drug design procedures.

Pharmacokinetics
Computational evaluation of the pharmacokinetic properties of the screened compounds using SwissADME showed high gastrointestinal absorption for the majority of the alkyl and alkenyl FA esters and chemically synthesized FA derivatives but low gastrointestinal absorption for the majority of the FA sugar esters (Table 6).

Toxicity Profile
Altogether, the compounds were predicted to have low toxicity, with LD 50 values being greater than 5000 mg/kg for the majority of them, and only a few alkyl, alkenyl, and fatty acid derivatives having a lethal dose lower than 1000 mg/kg. The prediction accuracy estimated by ProTox-II was around 70% for the vast majority of the compounds. Almost all appeared likely to be immunotoxic, but hepatotoxicity, carcinogenicity, mutagenicity, and cytotoxicity were predicted to be of no concern (Table 6).

Summary of the ADMET Properties of the Most Promising Compounds
Regarding the compounds with the best binding scores, namely FOS ferulate, FA rutinoside, galactobiose ferulate, raffinose ferulate, FA trimer, FOS ferulate 1, compound e27, xylobiose ferulate, arabinobiose ferulate, compound 4n, β-sitosteryl ferulate, compounds F3 and e28, and D-sucrose ferulate, the pharmacokinetic properties were not optimal. Among the best hits, only e27 has molecular weight, rotatable bond, polar surface area, logP, and logS values that indicate good oral bioavailability. The chemically synthesized derivatives, excluding 4n, which showed poor ADMET properties, are the only hits for which an adequate bioavailability score, good balance between water solubility and lipophilicity, and compliance with the Lipinski s rule were predicted. On the other hand, the enzymatically synthesized derivatives, comprising the sugar esters of FA, exhibit higher solubility and significantly lower acute toxicity, portrayed in the lethal dose values. Only compounds e27, e28, and F3 were predicted to have high gastrointestinal absorption while all of them appear to be unable to cross the BBB. Regarding selectivity, very few are assessed as cytochrome inhibitors, whereas none are a P-glycoprotein substrate.   6 : Prediction for the inactivity (green) or activity (red) of the compound in the respective toxicity category, as given by ProToxII; 7 : Probability of this prediction, given by ProToxII. An encouraging indicator is that most of the compounds are not P-glycoprotein or cytochrome inhibitors, meaning that they can be selective M pro targets, while they also seem to have acceptable skin permeability, taking into consideration that more negative values of logK p are interpreted as lower skin permeability [115]. As far as blood brain barrier (BBB) permeability is concerned, most compounds are not estimated to be able to permeate the BBB. This means that unwanted side-effects in the central nervous system are avoided. Nevertheless, there have been indications that SARS-CoV-2 does cross the barrier, potentially causing brain damage [118], so it is not entirely clear whether BBB permeability is an undesired property.
Based on all the above, FOS ferulate 1 would appear to have the highest balance between lipophilicity and solubility among the enzymatically synthesized derivatives, and an adequate druglikeness score, slight deviation from the oral bioavailability criteria (compared to the rest of the compounds), and a relatively high lethal dose. However, due to their higher binding scores, exceptional solubility, and the fact that they can be more easily synthesized, FA rutinoside and raffinose ferulate are also highlighted as promising hits. Regarding the chemically synthesized derivatives, compound e27 appears to have an edge over the rest of the compounds in all the categories, with the only concern being its relatively low lethal dose value. In any case, despite being a useful indication, these predictions are not definite and cannot exclude the compounds from further evaluation for their potential utilization as drugs or nutraceuticals.

MD Simulation of Selected Derivatives
MD simulation was performed to more thoroughly evaluate the stability of the complexes of three ligands that emerged as particularly promising from the molecular docking simulation (raffinose ferulate, FA rutinoside, and compound e27) with SARS-CoV-2 M pro while the complex of the protease with the native inhibitor N3 was also analyzed to provide a reference for the evaluation of the results. The complexes appeared to be particularly stable (Figure 12a), with overall low values and subtle fluctuation of the RMSD (1-2 Å) for FA rutinoside and compound e27, with the results for the latter indicating a complex that is more stable than that of N3 for the majority of the simulation time. Raffinose ferulate exhibited slightly higher RMSD values, which exceeded 2 Å after 6 ns. Although its RMSD values appear to have an ascending tendency, the simulation was run for 30 ns to confirm that equilibrium indeed occurred after 10 ns (results not shown). In any case, the RMSD values of the complex of raffinose ferulate are still below 3 Å, which is considered low enough to characterize the complex as stable [119]. The radius of gyration is a useful measure of the rigidity of the protein-l complexes. All the complexes appeared to have comparable results (Figure 12b), w can be regarded as indicators of compactness and stability [120,121]. Compoun showed the tightest complex with M pro . FA rutinoside had very similar RoG compa The radius of gyration is a useful measure of the rigidity of the protein-ligand complexes. All the complexes appeared to have comparable results (Figure 12b), which can be regarded as indicators of compactness and stability [120,121]. Compound e27 showed the tightest complex with M pro . FA rutinoside had very similar RoG compared to N3, and raffinose ferulate exhibited a greater variation in its RoG values, which, however, still compose a low RoG profile. Analysis of the RMSF values allows the acquisition of an image of which residues are more flexible. As seen in Figure 12c, the majority of residues have low RMSF values below 1 [120,122] while the terminal residues of the protease exhibit a much larger deviation and are the only ones to exceed 2.5 , with the exception of Arg222, Tyr154, Gln273, and Asp155 for compound e27. Therefore, the structure of the protein appears to be quite stable.
As far as protein-ligand interactions are concerned, the total hydrogen bonds formed between the two as a function of time are presented in Figure 12d and the number of total interactions between the ligands and each M pro residue (when occurring) are shown in Figure 13. The reference inhibitor N3 appears to be more stable in terms of its hydrogen bonds with the protease, which exhibit only slight deviations from the value of 3. Forming fewer hydrogen bonds but also showing good stability, compound e27 appears to have one hydrogen bond contact with the protease for most of the simulation time. The enzymatic FA derivatives FA rutinoside and raffinose ferulate form significantly more hydrogen bonds. In the case of FA rutinoside, the bonds appear to be more stable and fluctuate around the value of 6 while raffinose ferulate exhibits a more sudden drop after 2 ns of simulation to then be quite stable between 2 and 4 hydrogen bonds. Based on the number and stability of the hydrogen bonds, FA rutinoside shows indications of stronger binding to the protease.
Regarding the overall interactions, as seen in Figure 13a, FA rutinoside interacts with numerous residues around the active site, including multiple contacts with oxyanion hole residues Gly143, Ser144, and Cys145. Moreover, Val42 appears to interact with the ligand. Although it is not calculated as a contacting residue in the molecular docking simulation, it is a neighboring residue to the catalytic histidine (His41). Considering that the contacts are calculated based on the distance of the residue atoms from the atoms of the ligands, it can be difficult to identify and separate individual contacts of the ligand with neighboring residues. However, results show the vicinity of the ligand to this area of the active site. FA rutinoside also interacts with residues located at the S1 (Thr25, Thr26, and Asn28, which neighbor Leu 27, which is given as a contacting residue), S2 (Ser46, Glu47, Leu50), and S4 subsites (His164, Met165, Leu167). Raffinose ferulate also interacts with oxyanion hole residues, including the catalytic cysteine, and the residues Glu47 and Asp48, which are located near the contacting residues Thr45 and Met49 at the S2 subsite ( Figure 13b). Moreover, it interacts with several residues of the S4 subsite, with the ones appearing to form the most interactions being Thr190 and Ala193 (neighbor of the contacting residue Gln192), followed by Glu166. Lastly, compound e27 interacts with similar residues covering all four subsites (Figure 13c). Although it does not appear to interact with any of the catalytic residues, it forms contacts with their neighboring residues Val42, Ser144, and Gly146. Another observation is that the majority of the contacts are with residues located at the S4 subsite (Glu166, Leu167, Thr169, Arg188, Thr190, Gln192, Ala193).
Overall, compared to inhibitor N3 (Figure 13d), it is evident that the ligands examined form more contacts with the protease residues. The fact that the confirmed protease inhibitor N3 exhibits fluctuations in the number of contacts with the various contacting residues over time indicates that such fluctuations are normal. Although the contacts established between N3 and the protease appear to be more stable compared to the rest of the ligands, compound e27 also exhibits a consistent interaction profile. The contacts of FA rutinoside are also quite stable and higher in number overall. Raffinose ferulate, confirming the tendency shown from the RMSD and H-bond diagrams, appears to have the least balance in terms of interactions among the three candidates. In general, the MD simulation leads to the conclusion that the complex of compound e27 is the most stable among the three ligands investigated. However, the results yielded for FA rutinoside are also very encouraging and do not deviate much. This fact, combined with the greater number of interactions with the protease, establishes FA rutinoside as a very interesting candidate for the inhibition of M pro .
Biomedicines 2022, 10, x FOR PEER REVIEW 38 of 45 Figure 13. Number of contacts of the ligands (a) FA rutinoside, (b) raffinose ferulate, (c) compound e27, and (d) N3 with M pro residues throughout the simulation time. The residues marked with an asterisk (*) are residues that appeared to have contact with the ligand in the MD simulation but were not given as contacting residues in the molecular docking simulation output. However, they are neighboring contacting residues. This can be due to the fact that since MD trajectory analysis uses the distance as a parameter to calculate contacts, sometimes atoms from neighboring residues can be within a distance small enough to be regarded as contacts.
Regarding the overall interactions, as seen in Figure 13a, FA rutinoside interacts with numerous residues around the active site, including multiple contacts with oxyanion hole residues Gly143, Ser144, and Cys145. Moreover, Val42 appears to interact with the ligand. Although it is not calculated as a contacting residue in the molecular docking simulation, it is a neighboring residue to the catalytic histidine (His41). Considering that the contacts are calculated based on the distance of the residue atoms from the atoms of the ligands, it can be difficult to identify and separate individual contacts of the ligand with neighboring residues. However, results show the vicinity of the ligand to this area of the active site. FA rutinoside also interacts with residues located at the S1′ (Thr25, Thr26, and Asn28, which neighbor Leu 27, which is given as a contacting residue), S2 (Ser46, Glu47, Leu50), and S4 subsites (His164, Met165, Leu167). Raffinose ferulate also interacts with oxyanion hole residues, including the catalytic cysteine, and the residues Glu47 and Asp48, which are located near the contacting residues Thr45 and Met49 at the S2 subsite (Figure 13b). Moreover, it interacts with several residues of the S4 subsite, with the ones appearing to form the most interactions being Thr190 and Ala193 (neighbor of the contacting residue Gln192), followed by Glu166. Lastly, compound e27 interacts with similar residues covering all four subsites (Figure 13c). Although it does not appear to interact with any of the catalytic residues, it forms contacts with their neighboring residues Val42, Ser144, and Figure 13. Number of contacts of the ligands (a) FA rutinoside, (b) raffinose ferulate, (c) compound e27, and (d) N3 with M pro residues throughout the simulation time. The residues marked with an asterisk (*) are residues that appeared to have contact with the ligand in the MD simulation but were not given as contacting residues in the molecular docking simulation output. However, they are neighboring contacting residues. This can be due to the fact that since MD trajectory analysis uses the distance as a parameter to calculate contacts, sometimes atoms from neighboring residues can be within a distance small enough to be regarded as contacts.

Conclusions
The difficulties that have arisen in combating the ongoing pandemic highlight the importance of utilizing any available tool for immunity boosting. In this context, phytochemicals emerge as valuable allies. FA is an abundant bioactive compound, whose derivatization has been investigated both through chemical and enzymatic routes, in search for improved properties. This work studied 54 FA derivatives, among which 14 exhibited similar or better in silico binding affinity to the SARS-CoV-2 M pro compared to confirmed inhibitors. Further computational evaluation of their ADMET properties indicated FA rutinoside, raffinose ferulate, and compound e27 as the most promising hits to be further examined through molecular dynamics simulation. Analysis of the MD trajectories identified FA rutinoside and compound e27 as promising candidates, representing enzymatically and chemically synthesized derivatives, respectively. Several derivatives exhibited good binding affinity to the main protease though, indicating that it would be worth investigating different administration routes in order to overcome the problem of low oral bioavailability, and potential structural modifications that could lead to even higher affinity and more favorable pharmacokinetic properties. Moreover, the favorable results for many enzymatically synthesized derivatives encourage the development and optimization of more sustainable enzymatic processes, which can also include the valorization of biomass in which FA can be found. Even though further in vitro and in vivo research is necessary to provide more reliable data regarding the efficacy of the compounds, this work suggests that FA derivatives are interesting candidates to be considered for their antiviral potential against the current public health concern SARS-CoV-2 and future viral threats.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biomedicines10081787/s1, Figure S1: Binding modes of the second and third clusters for alkyl and alkenyl FA esters at the active site of M pro ; Figure S2: Binding modes of the second and third clusters for fatty acid, polyol, tocopherol, and sterol FA derivatives at the active site of M pro ; Figure S3: Binding modes of the second and third clusters for monosaccharide esters of FA at the active site of M pro ; Figure S4: Binding modes of the second and third clusters for disaccharide esters of FA at the active site of M pro ; Figure