Fourteen immunomodulatory alkaloids and two prenylated phenylpropanoids with dual therapeutic approach for COVID-19: molecular docking and dynamics studies

Abstract The pandemic outbreak of COVID-19 caused by the new severe acute respiratory syndrome coronavirus (SARS-CoV-2) is a global health burden. To date, there is no highly effective antiviral therapy to eradicate the virus; as a result, researchers are racing to introduce new potential therapeutic agents. Alternatively, traditional immunity boosters and symptomatic treatment based on natural bioactive compounds are also an option. The 3-chymotrypsin-like protease (3CLpro) crystal structure, the main proteolytic enzyme of SARS-CoV-2, has been unraveled, allowing the development of effective protease inhibitors via in silico and biological studies. In COVID-19 infected patients, the loss of lung function, and mortality are reported to be linked to several inflammatory mediators and cytokines. In this context, the approach of introducing immunomodulatory agents may be considered a dual lifesaving strategy in combination with antiviral drugs. This study aims to provide immunomodulatory natural products exhibiting potential protease inhibitory activities. Selected groups of alkaloids of different classes and two prenylated phenylpropanoids from the Brazilian green propolis were in silico screened for their ability to inhibit COVID-19 3CLpro protease. Results showed that compounds exhibited binding energy scores with values ranging from −6.96 to −3.70 compared to the reference synthetic protease inhibitor O6K with a binding energy score of −7.57. O6K binding energy was found comparable with lead phytochemicals in our study, while their toxicity and drug-likeness criteria are better than that of O6K. The activities of these molecules are mainly ascribed to their ability to form hydrogen bonding with 3CLpro crucial amino acid residues of the catalytic site. In addition, the molecular dynamics simulations further showed that some of these compounds formed stable complexes as evidenced by the occupancy fraction measurements. The study suggested that the major immunomodulators 3β, 20α-diacetamido-5α-pregnane, (20S)-(benzamido)-3β-(N,N-dimethyamino)-pregnane, and baccharin are 3CLpro inhibitors. Biological screenings of these phytochemicals will be valuable to experimentally validate and consolidate the results of this study before a rigid conclusion is reached, which may pave the way for the development of efficient modulatory bioactive compounds with dual bioactions in COVID-19 intervention. Communicated by Ramaswamy H. Sarma Graphical abstract was designed using BioRender (license #2702-5381)

Graphical abstract was designed using BioRender (license #2702-5381) Introduction SARS-CoV-2 is a new coronavirus, which was first identified in Wuhan and China in December 2019 as the causative agent of COVID-19.The infection caused by SARS-CoV-2 can lead to severe respiratory syndrome, which can lead to death (Kupferschmidt & Cohen, 2020) and the World Health Organization (WHO) and the Public Health Emergencies of International Concern declared COVID-19 as a pandemic (Kumar et al., 2023).This virus is a relatively large microorganism with a single-stranded RNA surrounded by an envelope consisting of a lipid bilayer and viral proteins.It belongs to beta coronaviruses, which mainly affects the lower respiratory system resulting in viral pneumonia.SARS-CoV-2 can also affect other organs such as the gastrointestinal system, kidney, heart, liver, and central nervous system.In severe cases, this can lead to multiple organ failure.(Su et al., 2016).
Several therapeutic interventions have been adopted for the management of COVID-19 including the use of antivirals, anti-inflammatory agents, and antithrombotics, among others (Murakami et al., 2023).However, there is no completely effective treatment for COVID-19.The virus can attack human cells using the spike protein (S protein), which binds to angiotensin-converting enzyme 2 (ACE2).Then, the S protein is cleaved into S1 and S2 subunits by a human cell-derived protease.The S1 binds to its receptor, ACE2 meanwhile the other fragment, S2 is cleaved by another serine protease enzyme (TMPRSS2), resulting in membrane fusion.The ACE2 and TMPRSS2 are the main cell entry enzymes involved in SARS-CoV-2 infection and spreading (Hoffmann et al., 2020).Inside the host cell, the viral genome is released and translated into viral polyproteins.These proteins are cleaved by a main viral proteinase (M pro ) into effector proteins.Then, a full-length negative-strand RNA template is synthesized to produce more viral genomic RNA (Liu et al., 2020).The M pro is often called 3-chymotrypsin-like protease (3CL pro ) and it is the main viral proteinase in the viral replication process.In addition, it is a promising target for the development and design of molecules inhibiting SARS and other coronaviruses infections (Anand et al., 2003).
Recent studies reported that there is a mild to severe cytokine storm in severe COVID-19 patients, which is suggested to be responsible for lung injury and pulmonary failure leading to death.Therefore, the management of cytokine storms and inflammation has become an important part of rescuing patients with severe conditions (C.Zhang et al., 2020a).In this regard, introducing immunomodulatory agents to modulate the elevated levels of cytokines can be potentially considered as a preventative approach to interfere with the disease progression in COVID-19 patients.
In the area of drug development, natural products represent more than 60% of the market pharmaceuticals (Molinari, 2009).Alkaloids are a major class of natural products that have been broadly used for the treatment of a wide variety of diseases (Amirkia & Heinrich, 2014;Cordell et al., 2001).They possess pharmacophores similar to antiviral drugs, which make them potential candidates as 3CL pro inhibitors (Fielding et al., 2020).In addition, several alkaloids, specifically the steroidal ones are previously reported as immunomodulatory agents (Gr€ uter et al., 2020;Iqbal et al., 2015;Kalmarzi et al., 2019;Peñaloza et al., 2023;Pathak & Khandelwal, 2009;Shen et al., 2020;Sunitha & Nagulu , 2018;Thawabteh et al., 2019).Some prenylated phenylpropanoids have been investigated due to their immunomodulatory and anti-inflammatory activities, which aroused our attention to screen atrepillin C and baccharin as potential antiviral candidates.An authentic green propolis extract obtained from Baccharis dracunculifolia, a plant native to Brazil, contains high content of artepillin C and baccharin which are the major and valuable bioactive principles found in the Brazilian green propolis (de Barros et al., 2007).Several studies have reported that green propolis extract and isolated major compounds can boost the body's defense mechanism and promote healthy body functions (Ferreira et al., 2021b).
The present study suggests repurposing fourteen immunomodulatory alkaloids and two prenylated phenylpropanoids with potential antiviral pharmacophores for the development of modulatory agents with a dual biological strategy for COVID-19 therapeutic intervention.This approach will not only reduce the cost and time of the drug screening process but also it may generate bioactive molecules with two dual actions; antiviral and immune system modulatory effects to prevent disease progression and decrease mortality rate.

3CL pro molecular docking study
The x-ray crystal structure of M pro (PDB ID 6y2f), the key protein involved in the virus life cycle, has been elucidated by Zhang et al. (2020b), which was reported with a covalently bound synthetic inhibitor O6K (Figure 1).This enzyme works by cutting the polyproteins translated from the viral RNA to yield functional viral proteins.Here, this study followed the same protocol regarding to selecting the active site of the protein for the grid preparation and the docking analysis.The grid box dimension and the grid center were set according to the centroid of the co-crystallized ligand structure (O6K), alpha ketoamide 13b (Zhang et al., 2020b).The compound O6K is an a-ketoamide inhibitor with a carbonyl warhead forming a C-S covalent bond with Cys145 (Mohapatra et al., 2023).The reported non-covalent and covalent docking energies of O6K were estimated to be À 7.118 kcal/mol and À 27.32 kcal/mol respectively (Mohapatra et al., 2023).

Ligand preparations
The three-dimensional structures of the selected compounds were downloaded from the Open chemistry database of PubChem in NIH https://pubchem.ncbi.nlm.nih.gov/.The ligand was prepared using Molecular Operating Environment (MOE) 2015.10 software (Chemical Computing Group, Montreal, Canada) (Corbeil et al., 2012).Partial charges were corrected by adjusting hydrogens and lone pairs as required prior to the energy minimization, then the potential energy was calculated.

Enzyme preparation
The crystal structure of the target protein (3CL pro enzyme) was downloaded from the protein data bank (www.rcsb.org) in PDB format with PDB code 6Y2E.The protein model was prepared using MOE software.The enzyme structure was three-dimensionally protonated after the addition of hydrogen atoms.The partial charge was calculated and MMFF94x force field was used.All ligands or water molecules were removed, and the site finder tool in MOE was used to select the binding sites based on information collected from the literature.

Docking and post-docking analysis
The ligand was docked in the prepared binding sites using the default settings of MOE.The default placement and refinement scoring functions are London dG and GBVI/WSA dG, respectively.The London dG scoring function determines the free energy of binding for a particular conformation position of the ligand.The GBVI/WSA dG Scoring function determines the free energy of binding for a particular conformation position of the ligand and is calculated using force fields of MMFF94x and AMBER99 (Corbeil et al., 2012;Love et al., 2022;Kurnikov et al. 2022).The post-docking Figure 1.The structure of the co-crystalized ligand O6K (alpha-ketoamide), the reference inhibitor, shows the important binding moieties involved in the interaction with 3CL pro crucial amino acid residues: Cys145, His41, and Glu166.
analysis provided the total energy required for binding and the docked poses were visualized in two-and three-dimensions using MOE.

Molecular dynamics simulations
The 3CL pro protein with PDB code 6Y2F (receptor) and the ligands (i) 3b, 20a-diacetamido-5a-pregnane, (ii) (20S)-(benzamido)-3b-(N,N-dimethyamino)-pregnane, (iii) baccharin, (iv) artepillin C, (v) capsaicin, (vi) piperine, (vii) tecostatine, (viii) 4-hydroxytecomine, (ix) tecomine, and (x) theophylline (ligands) were selected for the MD simulations.The co-crystalized ligand O6K was selected as a reference molecule in our experiments, as reported by (Zhang et al., 2020b).A system in the absence of any ligand (apoprotein) was also considered.The initial configuration of each system (receptorligand) was based on the outcomes given by the Molecular Docking simulations previously carried out.The systems were set up through the free online server CHARMM-GUI (http:// www.charmm-gui.org)(Jo et al., 2008).The parameterization of each ligand was performed with the Small Molecule Library (CSML) incorporated in CHARMM-GUI.In this way, the force field used for each ligand was CHARMM general force field (CGenFF) (Vanommeslaeghe & MacKerell, 2012).The force field for the receptor was CHARMM36m (Huang et al., 2017).All systems were solvated with the TIP3P water model and a monovalent salt (NaCl) concentration of 0.15 M. Additional ions were added in order to neutralize the systems.The GROMACS package version 2021.1 was employed for the MD runs (Abraham et al., 2015).The Particle Mesh Ewald (PME) method was used to treat long-range electrostatic interactions.Non-covalent interactions were set with a cutoff of 1.2 nm.The integration step employed was 2 fs.With the steepest descent algorithm, the energy minimization was performed until the systems reach a maximum force �1000 kJ mol À 1 nm À 1 .As usual, an equilibration phase was done in 310 K during 125 ps and imposing positional restraints to the backbone (400 kJ mol À 1 nm À 1 ).During this simulation, an NVT ensemble was employed and the Nose-Hoover algorithm was set with a coupling constant of s T ¼ 1.0 ps.The resulting equilibrated systems were used for the production phase maintaining the pressure at 1 bar through the Parrinello À Rahman barostat (with a coupling constant of s P ¼ 5.0 ps and compressibility of 4.5 � 10 À 5 bar).Each production run was simulated to reach 200 ns.

Measures for MD trajectories
By using the GROMACS package, two types of RMSD employing the rms tool: (1) were collected where the variations of the atomic coordinates of the protein were assessed and (2) where the motions of the ligands relative to the position of the protein were measured (labeled as 'RMSD lig ').For all cases, the reference structure was the initial configuration of the Docking outcomes.
To measure the frequency of interaction of each ligand with the protein surface, the occupancy fraction (OF) was computed with the program select of GROMACS with a cutoff ¼ 0.3 nm.When the OF is close to 1.0, it means that the interaction of the ligand with a specific region of the protein is occurring about 100% of the time in the trajectory.The occupied fraction (OF) was considered significant at a value � 0.5 according to (Barbera et al., 2018).For these estimates, a special emphasis on the pocket-targeted regions.In addition, the number of contacts per residue that occur between protein and the ligand over time has also been quantified.
From the MD trajectories, the binding energy was calculated from the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) using the g_mmpbsa program (Kumari et al., 2014).The total binding free energy (DG) for each case was calculated from the following equation: where G complex is the average free energy of the complex, and G 3CLpro and G ligand are the average energy values of the 3CL pro and ligand.The free energy of the protein complex and free protein is decomposed as follows:

Results and discussion
The treatment protocol for COVID-19 could be either therapeutic or preventative.Therapeutic intervention includes using direct-acting antivirals, host target agents, and immune-regulation therapy; meanwhile, preventative intervention employs vaccine development.The currently developed vaccines cannot present complete prevention from the virus due to continuous viral mutation.Part of the genetic code of the virus is changing during the viral translation process causing mutation.Therefore, other strategies targeting the virus via drug repurposing of novel biological pathways and modes of action could shorten the timeline for antiviral drug discovery and development of emergent pandemics.
According to some previous studies on 3CL pro enzyme and its crystal structure (Figure 2), the peptidomimetic alpha keto-amides are previously reported as potential 3CL pro inhibitors (Zhang et al., 2020b, Barbera et al., 2018).Through the nucleophilic attack of the catalytic Cys145 onto the a-keto group of an inhibitor, a thiohemiketal is formed in a reversible reaction.The oxyanion (or hydroxyl) group of this thiohemiketal is stabilized by a hydrogen bond from His41, whereas the amide oxygen of the inhibitor accepts a hydrogen bond from the main-chain amides of Gly143, Cys145, and partly Ser144, which form the canonical 'oxyanion hole' of the cysteine protease.One of the advantages of the a-ketoamides is that their warhead can interact with the catalytic center of the target proteases through two hydrogen bonding interactions rather than only one as with other warheads such as aldehydes or Michael acceptors (Zhang et al., 2020b).While in the Bi-cycloproline-containing 3CL pro inhibitors, the covalent bond is formed by Cys145 and the warhead aldehyde, which is essential for the antiviral activity as shown in (Figure 3) (Qiao et al., 2021).In addition, steroidal compounds are reported as NSP15 endoribonuclease inhibitors (Matsuyama et al., 2020), therefore they play a role in blocking the requisite viral entry process (Baglivo et al., 2020;Yao et al., 2020).Furthermore, they possess reported immunomodulatory activities (Horrigan et al., 2006;Jantan et al., 2015).
Collectively, the molecular fragments shown in (Figure 4) are favorable for SARS-CoV-2 3CL pro inhibitory activity, which is an attractive drug target against coronaviruses due to its vital role in processing the polyproteins, which are translated from the viral RNA (Amin et al., 2021).
The objective of this study is to identify potential antiviral natural compounds.Our research primarily explores alkaloidal compounds due to their significant structural features that enable them to bind with viral enzymes.To achieve this goal, we examined fourteen alkaloids with different chemical structures and classes (Figure 5), such as berberine (Abd El-Salam M. et al., 2015) , a benzylisoquinoline alkaloid piperine (a piperdine alkaloid), caffeine, theophylline, and theobromine (purine alkaloids), solanidine, capsaicin (a non-heterocyclic alkaloid), 3b,20a-diacetamido-5a-pregnane and (20S)-  (benzamido)-3b-(N,N-dimethylamino)-pregnane (steroidal alkaloids), and (boschniakine, 4-hydroxytecomine, N-normethylskytanthine, tecomine, and tecostanine), which are pyrindane (monoterpene) alkaloids.The selection of this group of alkaloids was based on two important major factors.The first is that they contain the potential pharmacophoric moieties for binding with 3CL pro .Second, they are reported as immunomodulators in several studies (Horrigan et al., 2006;Jantan et al., 2015), where they have the ability to inhibit the proinflammatory mediators responsible for provoking the cytokine storm in the respiratory system, which could worsen the condition of infected patients leading to death.Of note, steroidal alkaloids exhibit additional merit as promising agents for preventing viral entry (Andrade et al., 2017).We also have selected two major phenylpropanoids in the Brazilian green propolis extract: artipillin C and baccharin due to their reported potent immunomodulatory activities (Magnavacca et al., 2022).A previous study showed that baccharin mitigated inflammation by modulating the production of cytokines and eicosanoids (Ferreira et al., 2021b).
Remdesivir, molnupiravir, boceprevir, and telaprevir are the major antiviral drugs approved for treating COVID-19 infection.The UK Medicines and Health Products Regulatory Agency (MHRA) approved molnupiravir for treating adults with mild to moderate COVID-19 in 2021.Thus, molnupiravir was the first orally administered anti-SARSCoV-2 agent worldwide.Studies revealed that the SARS-CoV-2 Omicron variant is highly susceptible to molnupiravir (Li et al., 2022).However, the potential severe action of molnupiravir of provoking mutation in mammalian cells has raised major concerns (Lei et al., 2022;Zhou et al., 2021).Therefore, the search for safer and more effective anti-SARSCoV-2 agents is  (2021).Based on the structure of boceprevir and another peptidomimetic compound, the HCV inhibitor telaprevir, Qiao et al. developed new bicycloproline-containing M pro inhibitors (Qiao et al., 2021).A study performed by Xia et al. also showed the broad-spectrum anti-SARSCoV-2 action of two rationally designed inhibitors based on the peptidomimetic compounds telaprevir, and boceprevir (Xia et al., 2021).Other peptidomimetics such as MDL-28170, MG-132, and talaprevir, were also shown to possess M pro binding affinity.The cyclopeptide RTD-1, having anti-COVID-19 activities, was reported to be safe to support its investigation for the treatment of SARSCoV-2 infection (Park et al., 2022).Accordingly, the bicycloproline-containing 3CL pro inhibitors design was the rationale used in this present study, which prompted searching for safe compounds exhibiting these SAR.The selected phytochemicals were virtually screened for their ability to bind with 3CL pro using MOE (Corbeil et al., 2012).The results are represented in Table 1 and Figure 6.All compounds showed promising binding energy values with 3CL pro .Furthermore, they exhibited hydrogen-bonding interaction with 3CL pro at the crucial amino acid residues, Cys145, His41, and Glu166.According to previous studies, Cys145, His41 are mainly responsible for the enzyme catalytic activity, meanwhile, Glu166 adopts the conformation of the substrate binding site (Zhang et al., 2020b).
The compounds 3b, 20a-diacetamido-5a-pregnane, baccharin, and atripillin C were the top three ranked compounds possessing good scores compared to the reference inhibitor O6K.The molecule 3b, 20a-diacetamido-5a-pregnane exhibited hydrogen bond interaction with His41 amino acid residue with a binding energy score of À 6.53.In the case of baccharin, two hydrogen bonds were formed with the amino acid residues Cys145 and Thr26, in addition to a hydrophobic interaction with Gln189 with a score of À 6.96.Atripillin C showed a hydrogen bond and ionic interaction with His41 amino acid residue with a score of À 6.41.
Propolis is a viscous resin produced by Apis mellifera bees and has immunomodulatory activities (Berretta et al., 2017).Baccharis dracunculifolia DC (Asteraceae) is the most important botanical source of the green propolis native to southeastern.Brazil.Artepillin C and baccharin are the major prenylated phenylpropanoids (Figure7) previously isolated from the Brazilian green propolis (Watanabe et al., 2021).The major bioactive phenolic acid artepillin C is a diprenyl-phydroxycinnamic acid derivative.Numerous biological activities have been described for this compound, such as antiinflammatory, antimicrobial, antioxidant, antitumor, and gastroprotective (Beserra et al., 2021).
Tecoma stans alkaloids exhibited promising results as 3CL pro inhibitors except for N-normethylskytanthine, which showed no binding interactions with the amino acid residues in the active site.Tecomine exhibited two hydrogen bond interactions with His41 with À 4.79 score.Boschniakine and tecostanine formed hydrogen bonding with His163 and Cys145 with a score values of À 4.39 and À 4.67, respectively.Three hydrogen bonds were detected with the amino acid residues Cys145 and Glu166 with À 4.77 score in the case of 4-hydroxytecomine.
Even though Tecoma alkaloids except N-normethylskytanthine showed moderate binding scores, they are considered potential 3CL pro inhibitors due to their interaction with the crucial amino acid residues His41 and Cys145, and Glu166.

Molecular dynamics of 3CL pro and ligands
MD simulations are considered to be a reliable approach to investigating the underlying dynamic consequences of protein-ligand complexes (Ahmed et al., 2018;Alamri et al., 2021;Liu et al., 2017).MD simulations were performed for the selected compounds in the complex with 3CL pro (Figure S1).Thereby, MD simulation was performed based on ten ligands and the main protease.The first measure assessed was the root mean square deviation (RMSD).The RMSD squares between the structures created during the molecular dynamic simulation in the dimension of time is a common standard to confirm the protein structural stability alone (Apo conformation of the main protease) and in the presence of the ligand.This parameter indicates how the coordinates of a molecule change over time.Therefore, the RMSD of Ca atoms of the protein backbone was monitored throughout the 200 ns simulation (Figure S2).MD simulation of the inhibitor protease system indicates that the protease structure remains very similar to the X-ray structures (RMSD ranged from �0.1 to �0.4 nm, Figure S2).Overall, the main protease remained more stable in the presence of the bound inhibitors i.e. 3b, 20a-diacetamido-5a-pregnane, baccharin, artipillin C, and capsaicin when compared to its abo conformation.Similarly, the same pattern was observed for the co-crystalized ligand O6K, as illustrated in the yellow curve (Figure S2).On the other hand, (20S)-(benzamido)-3b-(N,Ndimethyamino)-pregnane, tecomine, and 4-hydroxytecomine where showing higher RMSD fluctuations due to important conformational changes are happening in these systems because of the presence of the ligand.
Calculations of the RMSD ligand (RMSD lig ) regarding the protein position also were computed (Figure S3).The first three compounds (i.e.3b, 20a-diacetamido-5a-pregnane, (20S)-(benzamido)-3b-(N,N-dimethyamino)-pregnane, and baccharin represented through the black, red, and blue curves respectively at the top of the plot) reach stability faster compared to the others.The reference molecule (i.e.O6K) was behaving similarly as shown in the yellow curve of Figure S3 even with less stability compared to the top three compounds (i.e.low RMSD lig value).After 70, 50, and 100 ns the compounds 3b, 20a-diacetamido-5a-pregnane (0.7 nm), (20S)-(benzamido)-3b-(N,N-dimethyamino)-pregnane (1.5 nm), and baccharin (1 nm) reach a configurational equilibrium, respectively.Note that the values in parentheses refer to the RMSD lig value achieved for each ligand.The first two ligands belong to steroidal alkaloids and the remaining third is a prenylated phenylpropanoid.Other organic molecules such as artepillin C, capsaicin, and tescostatine also reach configurational equilibrium, but they were far from the pocket of interest.They are fluctuating too much as shown in Figure S3.It is important to highlight that Tecoma stans alkaloids (pyridine (monoterpene)) and purine alkaloids contain some potential key structural units for SARS-CoV-2 inhibition.Nevertheless, these molecules have not shown meaningful stability in the MD study.We speculate that this may be because the length of their structure does not fit in the binding pocket or non-covalent interactions are not optimized enough to keep the integrity of the complex.
The occupied fraction (OF) for the ten selected ligands interacting on the surface of 3CLpro was calculated.A OF greater than 0.5 is showing a meaningful association between the ligand molecule and a specific protein residue.Furthermore, it also means that the interaction between two specific molecules is happening half of the time overall trajectory.This study emphasizes the binding target pocket as shown in Figures 8 and S1 (the highlighted regions in cyan and yellow, respectively).Based on these criteria (i.e.OF > 0.5 and proximity to the binding pocket) the most relevant ligands were mainly the three molecules 3b, 20a-diacetamido-5a-pregnane, (20S)-(benzamido)-3b-(N,N-dimethyamino)-pregnane, and baccharin compared to the reference compound O6K.The other remaining compounds did not show significant OF values or interacted very far from the pocket of interest (Figure 8).
The three molecules with the best OF values are in good agreement with our obtained docking results.This study highlighted that the pose predicted for the docking assays was relatively well preserved specially for the MD simulation of 3b, 20a-diacetamido-5a-pregnane (mean RMSD lig ¼ (0.6 ± 0.2) nm).The remaining two molecules with good OF (i.e.(20S)-(benzamido)-3b-(N,N-dimethyamino)-pregnane, and baccharin) showed different configurations with respect to the docking prediction but stayed close to the binding pocket with important stability (check also the red and blue RMSD lig -curves, respectively, at the top plot in Figure S3).These three compounds interact with His41, Cys145, Glu166, Gln189, and Gln192 amino acid residues in the targeted binding pocket, which are crucial for binding with 3CL pro enzyme.Among these residues, Cys145 and His41 are relevant because they are considered catalytic dyads which is an important feature for the inhibition of protease activity (Ferreira et al., 2021a).Thereby, the three compounds 3b, 20a-diacetamido-5a-pregnane, (20S)-(benzamido)-3b-(N,N-dimethyamino)-pregnane, and baccharin could be potential therapeutic/drug candidates for enzyme viral inhibition.
The number of contacts between the first three compounds and the reference compound O6K with better OF values and the enzyme are ranging from 10 to 15 contacts.This number of contacts is maintained stable along the trajectory (Figure 9).Moreover, molecules such as artepillin C, capsaicin, and tecostanine present relative stability at the level of the number of contacts in potential regions that could be pocket candidates for enzymatic inhibition despite low OF (<0.5).It must be confirmed in the residue regions 248-260 and 290-308 for capsaicin, tecostanine, tecomine; and 270-290 for artepillin C shown in Figure 8.The compounds (i.e. piperine,tecomine,and theophylline) show fluctuations regarding the number of contacts, and they also explore different regions of the protein surface, which indicates a low affinity for the pocket of interest (Figures 8 and S1).
It is important to highlight that the most unstable molecules in terms of the interaction with the target binding pocket are also the smallest (artepillin C ¼ 22, capsaicin ¼ 22, piperine ¼ 21, tecostatine ¼ 13, 4-hydroxytecomine ¼ 14, tecomine ¼ 13, and theophylline ¼13) compared to 3b, 20a-diacetamido-5a-pregnane (¼ 30), (20S)-(benzamido)-3b-(N,N-dimethyamino)-pregnane (¼ 35), and baccharin (¼ 27), and the reference O6k (¼ 43).The numbers in parentheses show the number of heavy atoms per compound.These results may be showing an important precedent and we hypothesize that is probable that a number of sites greater than 27 are necessary for stable interaction with 3CLpro (at least for the ligand molecules explored here).Cys145, Glu166 -À 7.57 ± 0.08 � Despite being prenylated phenylpropanoids, these compounds have shown promising in silico results similar to the main alkaloids we have investigated.They are major constituents of Brazilian green propolis extract, which is known for its diverse potential biological activities, including antiviral effects against herpes simplex virus type 1 (HSV-1) and others (Coelho et al., 2015).
However, the stability of these couplings during molecular dynamics is strongly dependent on non-covalent interactions such as H-bonds, salt bridges, or van der Waals forces (vdW), and protein conformational changes.Therefore, the ligand number of sites hypothesis proposed here should be considered in better detail.An evaluation of the van der Waals interactions and hydrophobic associations could give a clue due to most of the residues in the target pocket being essentially hydrophobic.Thereby, a high number of sites in the ligand could have a direct relationship with the coupling stability since as far as the surface molecule is higher the vdW forces also could be increasing (Israelachvili, 2011).To test this hypothesis, it was necessary to calculate the binding free energy of protein-ligand complexes in order to assess the contribution given for short-or long-range interactions.

MM-PBSA free binding energy calculations
The MM-PBSA is a popular method for predicting free energy of binding due to its good accuracy compared to most other scoring functions of molecular docking methods (Wang et al., 2019).MM-PBSA based binding free energy of all protein-ligands complexes was calculated for 200 ns trajectory.The binding free energy was determined using MMPBSA method.Furthermore, through this approach, it is possible to decompose in different energetic terms the electrostatic, polar solvation, van der Waals, and SASA contributions that drive the interaction between the 3CLpro and the ligands (Table 2).The reference compound O6K showed free binding energy of À 65.79 kJ/mol whereas 3b, 20a-diacetamido-5apregnane, (20S)(benzamido)-3b-(N,N-dimethyamino)-pregnane and baccharin showed binding energy of À 67.81, À 90.07, and À 80.93 kJ/mol, respectively.These compounds are showing higher binding affinity for 3CL pro compared to O6K.
Regarding the three first compounds and O6K, it is observed that the forces are contributing significantly to the association with the viral protein compared to the rest of the ligands (< À 127 kJ/mol).In parallel, the SASA energies of the first three ligands and O6K showed more negative values compared to the remaining 8 ligands (< À 16 kJ/mol).This could support the previously raised hypothesis of this study.It is worth mentioning that the polar solvation energies are also the most positive among these four ligands, specifically in the case of compound baccharin and the reference inhibitor O6K.Other compounds which stood out with negative total binding free energies were Artepillin C (À 42.75 kJ/mol) and Capsaicin (À 62.53 kJ/mol).In the case of these compounds, it was clear that vdW forces, despite being far from the target pocket, governed the net attraction.The remaining compounds (piperine, tecostatine, and theophyline) exhibited low affinity values.
Since the outbreak of COVID-19, many computational studies have been conducted to unveil new potential 3CL pro inhibitors from diverse sources including FDA-approved drugs, natural products, synthetic small molecules, and peptides.For instance, some in silico screening studies have identified novel inhibitors from marine products (Gentile et al., 2020), flavonoids (Batool et al., 2022), protease inhibitors (Keretsu et al., 2020), and commercial chemical libraries (Ibrahim et al., 2021).To the best of our knowledge, there is no previous study covering the rationale of the dual action of the natural products was reported.The current study was based on structure similarity studies between the selected fourteen alkaloids and the two prenylated phenylpropanoids with previously reported reference inhibitors as mentioned above.
Virtual screening through molecular docking has some limitations including variability in the predicted scores (Koes et al., 2013).To overcome these limitations, this study repeated the molecular docking experiments to verify and validate these obtained results.Moreover, the present study analyzed protein-ligand interactions through MD simulations and free energy calculations using MM-PBSA method and applying different important parameters such as RMSD, OF, and the number of contacts per residue.Conceivably, our insilico study could be an adjunct to, not a substitute for, experimental validation of inhibitors for SARS-CoV-2 3CL pro .

Conclusion
Our structure-based virtual screening and molecular dynamic study conducted on selected immunomodulatory phytochemicals reveals promising 3CLpro inhibitors.Both steroidal alkaloids and prenylated phenylpropanoids show potential, particularly the three phytochemicals 3b, 20a-diacetamido-5a-pregnane, (20S)-(benzamido)-3b-(N,N-dimethyamino)-pregnane, and baccharin, when compared to the reference inhibitor O6K.While these results require further validation through direct inhibitory assays using SARS-CoV-2 3CLpro, as well as SARS-CoV-2 in vitro cell culture assays and other preclinical antiviral screening models, our study speculates that these compounds may inhibit crucial steps in the viral replication cycle.Additionally, the immunomodulatory effects of these compounds may reduce the cytokine storm induced by SARS-CoV-2 infection, making them unique and having additional merit.This dual mode of action may pave the way

Figure 2 .
Figure 2. The crystal structure of 3CL pro .(a) Cartoon representation of the structural domains of 3CL pro of SARS-CoV-2 (PDB code 6Y2E).Domain I (residues 10-96) is shown in yellow, domain II (residues 102-180) in green, and domain III (residues 200-303) in pink.The active site is located at the interface between domains I and II.(b) The active site of 3CL pro showing the catalytic residue Cys145, which is part of domain II, is 3.6 Å from His41 of the catalytic dyad, which is part of domain I (Ferreira & Rabeh, 2020).

Figure 3 .
Figure 3. (a) Schematic diagram of bicycloproline-containing 3CL pro inhibitors design derived from the two FDA-approved antiviral drugs: (b) Chemical structures of boceprevir or telaprevir.These derivatives also contain favorable fragments for SARS-CoV-2 inhibition.The carbon of the warhead aldehyde interacts with the sulfur atom of catalytic residue Cys145 to form a 1.8-Å covalent bond.The oxygen of the aldehyde forms two hydrogen bonds with the main-chain amides of Cys145 and Gly143 (forming the 'oxyanion hole').The P1 lactam ring inserts deeply into the S1pocket.The oxygen and nitrogen of lactam form two hydrogen bonds with the side chain of His163 (2.8 Å) and the main chain of Phe140 (3.3 Å), respectively.The main-chain amide of P1 makes a 2.9-Å hydrogen bond with the backbone O of His164.Figure was adapted according to (Qiao et al., 2021).

Figure 4 .
Figure 4.Chemical structures of the molecular fragments of the screened phytochemical compounds (pharmacophoric moieties), which are favorable for SARS-CoV-2 3CL pro inhibitory activity.

Figure 5 .
Figure 5.Chemical structures of the selected alkaloids as potential promiscuous therapeutic agents against SARS-CoV-2.The outlines are showing the important pharmacophoric moieties (key structural units), which are proposed to be responsible for their antiviral activity.

Figure 6 .
Figure 6.The 2D molecular model of the selected compounds binding with 3CL pro obtained by MOE software showing the amino acids involved in the interaction.These compounds were docked into the active site of 3CL pro with PDB code 6Y2E.

Figure 7 .
Figure 7.Chemical structures of artepillin C and baccharin, the prenylated phenylpropanoides isolated from Brazilian green propolis as potential promiscuous therapeutic agents against SARS-CoV-2 showing the pharmacophoric moieties (key structural units) are considered to be responsible for their anti-viral activity.

Figure 8 .
Figure 8. Occupied fraction (OF) of the different phytochemical ligands interacting at different protein regions in comparison with O6K.The cyan lines in the background indicate the target pocket.

Figure 9 .
Figure 9. Number of contacts per protein residue with the different phytochemical ligands is compared to the reference synthetic compound O6K, as a function of time.

Table 1 .
Virtual screening results of the selected compounds against 3CL pro .

Table 2 .
Binding free energy (kJ/mol) calculations for the selected ligands and the reference ligand calculated by MM/PBSA.