Molecular modelling and structure-activity relationship of a natural derivative of o-hydroxybenzoate as a potent inhibitor of dual NSP3 and NSP12 of SARS-CoV-2: in silico study

Abstract The nsp3 macrodomain and nsp12 (RdRp) enzymes are strongly implicated in the virulent regulation of the host immune response and viral replication of SARS-CoV-2, making them plausible therapeutic targets for mitigating infectivity. Remdesivir remains the only FDA-approved small-molecule inhibitor of the nsp12 in clinical conditions while none has been approved yet for the nsp3 macrodomain. In this study, 69,067 natural compounds from the IBScreen database were screened for efficacious potentials with mechanistic multitarget-directed inhibitory pharmacology against the dual targets using in silico approaches. Standard and extra precision (SP and XP) Maestro glide docking analyses were employed to evaluate their inhibitory interactions against the enzymes. Four compounds, STOCK1N-45901, 03804, 83408, 08377 consistently showed high XP scores against the respective targets and interacted strongly with pharmacologically essential amino acid and RNA residues, in better terms than the standard, co-crystallized inhibitors, GS-441524 and remdesivir. Further assessments through the predictions of ADMET and mutagenicity distinguished STOCK1N-45901, a natural derivative of o-hydroxybenzoate as the most promising candidate. The ligand maintained a good conformational and thermodynamic stability in complex with the enzymes throughout the trajectories of 100 ns molecular dynamics, indicated by RMSD, RMSF and radius of gyration plots. Its binding free energy, MM-GBSA was recorded as –54.24 and –31.77 kcal/mol against the respective enzyme, while its structure-activity relationships confer high probabilities as active antiviral, anti-inflammatory, antiinfection, antitussive and peroxidase inhibitor. The IBScreen database natural product, STOCK1N-45901 (2,3,4,5,6-pentahydroxyhexyl o-hydroxybenzoate) is thus recommended as a potent inhibitor of dual nsp3 and nsp12 of SARS-CoV-2 for further study. Communicated by Ramaswamy H. Sarma


Introduction
The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a b-CoV member of the family coronaviridae, newly emerged in late 2019. Its infection was later termed coronavirus 2019 . Since its emergence, the world has been experiencing an unprecedented destabilization in all aspects of life, whilst the search for effective control strategy lasts. The effective therapeutic interventions for the prophylaxis and treatment of single-stranded, positive, ribonucleic acid (RNA)-genomic viral infections including the COVID-19 involve the targeting of viral entry and intracellular processes such as mutation, replication and transcription V'kovski et al., 2021). The replication and transcription of the SARS-CoV-2 RNA genome occur mainly through the RNAdependent RNA polymerase (RdRp) enzyme and the non-structural proteins (nsps) (Elkarhat et al., 2020;Kokic et al., 2021). Upon entry mediated by the angiotensin-converting enzyme 2 (ACE 2) and transmembrane protease serine 2 (TMPRSS2), the viral genome functions directly as a messenger RNA, translating the open reading frames 1a and 1 b (ORF 1a and ORF 1 b) polyproteins. The progressive cleavage of the two polyproteins produces sixteen nsps (nsp 1-16) ( Figure 1) which are strongly implicated in the assembly of virion for various virulence activities including the disruption of host immune response, genome replication and transcription (Brant et al., 2021;Kim et al., 2021;Ni et al., 2021).
The nsp3 macrodomain (also known as mac-1/EC 3.4.22.46) is an integral protein, encoding a papain-like protease (PL pro ) essential for the viral proliferation and polyprotein processes (protease activities) for the formation of nsp 1-3. Other nsps (4-16) are hypothetically processed by the chemokine-like protease 3 (3CL pro ), present in nsp5 (EC 3.4.22.69 (Ni et al., 2021;Osipiuk et al., 2021). Besides the PL pro domain of the nsp3, other catalytic domains include the ubiquitin-like (Ubl), adenosine diphosphate (ADP) ribose phosphatase (ADRP), SARS-unique domain and the second Ubl co-joined with the PL pro along the N-terminal, while the nucleic acid binding domain (NAB), transmembrane domain (TM) and zinc finger motif (ZnF) constitute the C-terminal (Figure 2(A)) (Armstrong et al., 2021). The nsp3 is the largest among the nsps and an essential component of the viral complex for replication and transcription assembled on the host membranes. The ADP ribose binding pocket of the protein consists of loops connecting the b3-a2, b5-a4 and b6-a5 helices lining the site for flexibility to accept various ligand inhibitors (Figure 2(B,C)) (Ni et al., 2021). The nsp3 is strongly implicated in membrane association and modification by coronaviruses including the SARS-CoV-2 such as the removal of posttranslational ADP ribosylation sites to regulate the host immune response. The posttranslational modification is usually in form of mono-adenosine diphosphate (ADP)-ribosylation (MARylation) or poly-ADP-ribosylation, catalysed by poly-ADP-ribose polymerase (PARP). The catalytic process is importantly associated with host immune response, impair of host defence and stress signalling. These make the protein crucial for viral existence, virulence activity and an important therapeutic target (Denison, 2008;Kim et al., 2021;Ni et al., 2021;Osipiuk et al., 2021).
The nsp12 is an RdRp enzyme for catalytic RNA replication through RNA synthesis. The nsp12-synthesised RNA enhance the catalytic integrity of other nsps for the production of genomic and sub-genomic RNAs (gRNA and sgRNA) outside the double membrane vesicles of the viral replication and transcription complex (RTC). The sgRNA in the infected cell thus, serves as the messenger RNA, essential for the production of accessory and structural proteins. The core RTC consists of the nsp12 pivotal for the viral existence, replication and transcription through RNA synthesis. In addition to the nsp12, the RTC core contains of three accessory subunits, consisting of one nsp7 and two nsp8. These are connected with nsp13(adenosine triphosphate-dependent 5 0 and 3 0 RNA helicases) which hypothetically promotes the RNA template while nsp9 modulates the cap synthesis (Figure 3(A)). The nsp3 additionally binds the nsp4 to form an intermediate of the RTC, thereby contributing to the cleavage of the large polypeptide product of ORF. Some antiviral agents including favipiravir, remdesivir and suramin have been demonstrated with inhibitory binding to the RTC domains ( Figure 3(B,C)) 1ab (Brant et al., 2021;Kim et al., 2021). As such, the binding of natural product inhibitors to the active sites of the dual targets, nsp3 and nsp12 represent effective multi-target therapeutic pathways for mitigating the viral virulence and polymerase activities for replication and transcription as shown in this study.
Although, the emergence of effective vaccines such as the Oxford-AstraZeneca, Pfizer-BioNTech and Moderna has aided a great relief globally from the pandemic of SARS-CoV-2 infections including the deadly Delta variant. However, some recent studies conducted across some countries including Qatar and UK have shown that their protection integrity wanes over time (Pouwels et al., 2021;Tang et al., 2021). Again, most of the currently available ones are encoded with modified-RNA (mRNA) of SARS-CoV-2 with mechanistic potentials for activating the human immune response, supporting the RdRp as a plausible therapeutic target. Interestingly, some have recorded >95% protective effects in vaccinated people, however, the development vaccine-escaping mutations in the viral S-protein also threatens their long term efficacy and dependability (Ibrahim et al., 2021b(Ibrahim et al., , 2021c. These make the search for effective anti-SARS-CoV-2 drugs remains crucial. Notably, remdesivir is an outstanding candidate, expressing the mechanistic inhibition of RdRp and remains the only successful FDA-approved nucleoside in clinical situations of COVID-19, acting through the target (Elkarhat et al., 2020). Being a phosphoramidate prodrug, its cellular metabolism produces remdesivir triphosphate (RTP). The RTP is utilized as a substrate by RdRp to incorporate remdesivir monophosphate (RMP) in place of adenosine monophosphate (AMP) for the RNA growth, which extends with three additional nucleotides. The binding of remdesivir induces the fourth nucleotide impaired by further RNA translocation barrier which prevents the entry of more nucleoside triphosphate and eventually stalls the RdRp (Kokic et al., 2021). Several bioactive molecules including poly(ADP-ribose) polymerase inhibitors have been reported with efficacious inactivation of the nsp3 macrodomain. However, most of these agents re macromolecules. Some recent molecular science studies identified small molecules, remdesivir active metabolite, GS-441524 and some b-carboline alkaloids as dual plasticity enhancers and inhibitors of the nsp3 macrodomain Brosey et al., 2021;Ni et al., 2021), creating opportunities for investigating other diverse small active molecules.
Multi-target pharmacology has become invaluable in modern medicinal science. The approach, "one-drug-multiple-target" has enhanced therapeutic breakthroughs to conquer some diseases with complex etiologies such as neurodegenerative disorders, cancer, microbial and viral infections (Talevi, 2015). Relevantly, several scientific investigations have resulted in the proposal of some potent inhibitors with multi-target mechanisms of action amenable for overcoming the drug resistance and comorbidities that trigger severity in COVID-19. These include silibinin, a natural product inhibitor targeting the dual inflammatory cytokine signalling of the host and viral replication process (Bosch-Barrera et al., 2020). Some antimicrobial agents including streptolydigin and rifabutin, interactively inhibit the nsp12 motifs and nsp7-nsp8-nsp12 complex structures (Elkarhat et al., 2020). Other natural products such as d-viriferin, myricitrin and hesperidin reportedly bind strongly to the viral main protease (M pro ), RdRp (nsp12) and host ACE2, indicating multi-target pharmacology against the SARS-CoV-2 (Joshi et al., 2020), and as repurposable drugs targeting the dual 3CL pro and PL pro (Rajpoot et al., 2021). Another documented research also supports the plausibility of dual targeting of the toll-like receptor 4 and ACE2 as an effective therapeutic pathway against the SARS-CoV-2 infectivity and pathogenesis (Kate Gadanec et al., 2021). The interesting outcomes of the studies inspire further exploration of potential inhibitors with multi-target pharmacology against some essential therapeutic targets of SARS-CoV-2 including the dual nsp3 and RdRp. Similarly, the databases containing safer and cheaper natural product prospects to fulfil the quest remain underexplored. The IBS INTERBIOSCREEN (IBScreen, https://www.ibscreen.com/) remains one of the world's leading providers of high quality natural/synthetic chemical libraries of for screening. Moreover, some recent relevant studies have scanned the database and identified promising SARS-CoV-2 inhibitors acting on various targets (Elekofehinti et al., 2021a(Elekofehinti et al., , 2021bThomas, 2021), as such found worthy of exploration.
Moreover, the application of computer-aided molecular modelling in modern drug design has aided the rapid and successful discovery of novel candidates with improved pharmacology amenable for wet-lab analyses. Molecular docking, quantum molecular mechanics and functional group structure-activity relationship (SAR) have been extensively applied to investigate approximated molecular mechanisms of disease targets in biomedical research (De Ruyck and Brysbaert, 2016;Harrold & Zavod, 2013). The application of the in silico techniques in some recent studies have also afforded the identification of some potent inhibitors of the SARS-CoV-2 main protease for experimental validation (Ibrahim et al., 2020(Ibrahim et al., , 2021a(Ibrahim et al., , 2021b(Ibrahim et al., , 2021c. Thus, this study is aimed at identifying potent inhibitors of the SARS-CoV-2 polymerase activities from natural product databases acting through the dual nsp3 and RdRp enzymes, using computational molecular modelling and in silico SAR study. The study extensively incorporates molecular docking, molecular dynamics (MD), in silico SAR study, predictions of absorption, distribution, metabolism, excretion, toxicology (ADMET) and mutagenicity.

Ligand retrieval and preparation
The chemical library of IBScreen (https://www.ibscreen.com/) was searched for screening compounds and 69,067 natural product compounds were retrieved in structure data file (SDF) format. The IBScreen database was selected because among the natural product databases, it is enriched with openly accessible and freely downloadable formats of lead-like small molecules. Some recent studies have scanned the database to identify promising SARS-CoV-2 inhibitors acting on various targets (Elekofehinti et al., 2021a(Elekofehinti et al., , 2021bSorokina & Steinbeck, 2020). The retrieved files include those of the references, remdesivir and its GS-441524 metabolite which were obtained from the PubChem database were prepared using the LigPrep module of Maestro 12.12 (LigPrep, Schrodinger, LCC, New York, NY, 2019). The protocols involved the assignment of protonation states, energy minimization, geometry and partial atomic charge estimation using the optimized potentials for liquid simulations 3e (OPLS-3e) force field (Harder et al., 2016). The Epik at a target was used in default settings to generate possible ionization states at pH 7.0 ± 2.0. Finally, the minimized and protonated ready-to-dock conformers were saved in LigPrep.out file.
Protein structure retrieval, preparation and active site cavity The RCSB protein data bank (PDB) (https://www.rcsb.org) was searched for protein structures of the SARS-CoV-2 nsp3 and nsp12 with respective resolution of 2.15 and 3.10 Ð, ideal for accurate representation of the targets for the study. The crystal structure of the nsp3 macrodomain in complex with GS-441524 (PDB 7BF6) and the electron microscopic structure of the elongating RdRp in complex with remdesivir (PDB 7B3B) were retrieved from RCSB in PDB formats. The structures were imported into the workspace of Maestro 12.2, preprocessed by assigning the bond orders, protonation states and filling of missing loops using the Protein preparation wizard (Madhavi Sastry et al., 2013). These were followed by optimization, refining and energy minimization through which the conversion of heavy atoms to relative mean standard deviation (RMSD) 0.3Ð was achieved using the OPLS-3e force field (Harder et al., 2016). The enclosing grid boxes within which ligands were expected to be actively docked were generated by selecting the centroid of the workspace ligand, GS-441524 and remdesivir. The x, y, z coordinates of the respective grid boxes for PDB 7BF6 and PDB 7B3B were obtained at default settings of the ligand docking length of 0-20 Å. The active site amino acid residues were identified and the pocket cavities, represented by the ready-to-dock receptor grid boxes were saved in gridbox.zip files.

Molecular docking
The inhibitory interactions of the 102,174 natural product protonated conformers generated from the LigPrep protocol were evaluated in silico against the respective targets in comparison with standard inhibitors, remdesivir and its metabolite, GS-441524 using ligand-receptor molecular docking. The GS-441524 and remdesivir co-crystallized in the respective structure of the nsp3 macrodomain (PDB 7BF6) and RdRp (PDB 7B3B) were used as controls. Both output files of the prepared ligands and grid box of each receptor file were uploaded into the Maestro 12.2 Glide docking panel for standard and extra precision (SP and XP) ligand-receptor docking analysis (Halgren et al., 2004). The docking panel were set at default to generate 5000 poses per ligand at the initial phases of the docking while best 400 poses were set to be produced for energy minimization using OPLS-3e. The distance-dielectric constant was set at 2.0 with a 100 maximum number of minimization steps. The SP module was first applied to screen a large number of the ligands while the top-ranked 100 ligands with good SP Gscores and binding poses were selected for XP docking. The XP module for the top-ranked ligands from the SP protocols provides an opportunity for more discriminate, robust torsional refinement and ligand sampling (Parmar et al., 2021). The precision, accuracy and reproducibility of the docking protocols were validated through redocking of the co-crystallized ligands, remdesivir and GS-441524 within the active pockets of the enzymes and computing their RMSD in comparison with the retrieved crystal complex from RCSB PDB. The RMSD values of 2.00 Å validates the protocols as ideally reliable (Castro-Alvarez et al., 2017;Ram ırez and Caballero, 2018). The docked ligands were ranked based on the algorithms of the SP and XP docking Gscores, and their interactions with pharmacologically essential amino acid residues in the active binding pockets of the enzymes analysed using the generated 2D binding poses. The natural product ligands with top-notch good scoring functions, binding poses and binding interactions, representing promising inhibitors relative to remdesivir and GS-441524 were mapped for further studies.

Prediction of pharmacokinetics, druggability and toxicological properties
The physicochemical, pharmacokinetics, pharmacodynamics and toxicological properties of the four top-ranked natural products from the docking analyses were predicted using the pkCSM-pharmacokinetics server by the input of SMILE files. This essentially screens the candidates with the most Figure 3. Structures of SARS-CoV-2 replication and transcription complex (RTC). (A) A composite structure of RTC from three PDB coordinates: PDB 7CXM architecture of nsp7, nsp8 X2, nsp12 and nsp13 X2 bound to RNA template and primer), PDB 6XEZ (the ADP・AlF3, bound in the nsp13 helicase active site and PDB 7CYQ (nsp9 associated with nsp12 and GDP in the active site of RNA capping). The RNA template pieces bound to nsp13 and nsp12 are not connected and would be pulled by the two enzymes in opposite directions as indicated by the yellow double arrowheads. (B) C Zoom-in views of RTC bound to inhibitors, favipiravir (PDB 7AAP), remdesivir (RMP) (PDB 7B3B), and suramin (PDB 7D4F). The nsp12 is shown in grey. The three inhibitors are in distinct colours, mimicking the phosphate backbone of RNA, two suramin molecules (cyan) compete for the RNA template and primer binding. remdesivir (blue) is already incorporated in the RNA primer strand at -3 position. Favipiravir RTD (magenta) occupies the incoming nucleotide position, but the phosphates are in a non-productive conformation. The active site residues are shown in pink-red sticks and Mg2þ ions are shown as green spheres. (Reprinted from (Brant et al., 2021), Licensed under CC-BY 4.0). (C) Binding interaction showing hydrogen bonding and p-p interaction of remdesivir with amino acid residues and RNA molecules at the active site along chain A of nsp12 (PDB 7B3B).
probable potentials for the delivery to targets for physiological expressions and removal from the system without acute toxicity. The module expresses a distance-based graph structural signature through accurate and precise calculations of atomic pharmacophore frequency counts, toxicophore fingerprints and general molecular properties for new candidates under probe for drugability (Pires et al., 2015). The predicted parameters include molecular weight (MW), number of hydrogen bond acceptor/donor (HBA and HBD) lipophilicity (Log P), number of rotatable bonds (for flexibility), solubility (Log S), topological polar surface area (TPSA). These and others quantitatively define the potentials of the selected candidates for ADMET. Their druggability was also determined using the parameters for blood-brain barrier (BBB), human gastrointestinal absorption, pglycoprotein (pg) substrate and the inhibitory potentials on the cytochrome P450 isoenzymes. The parameters were benchmarked with the BDDC rules of 5 (RO5) and druggability, and Veber's influence of molecular properties on oral bioavailability (Benet et al., 2016;Veber et al., 2002). The natural product with promising pharmacokinetics, druggability and toxicological profiles was identified for further studies.

Mutagenicity
The toxicological indices of the finally selected natural product, STOCK1N-45901 were further assessed through the prediction of mutagenicity using the read-across SAR-inclined VEGA ToxRead 0.23Beta java-based software. The fragment and molecular similarity search from the series of known compounds were modelled to predict the mutagenicity of the selected molecule, inputting the string of its SMILE format. The accuracy and reproducibility of the obtainable results have been validated with the Ames in vitro experiments for applicability (Gini et al., 2014).

Molecular dynamics
The best ranked with the most promising ADMET profile, STOCK1N-4590 in complexes with nsp3 macrodomain and RdRp (PDB 7BF6 and PDB 7B3B) was further studied for assessing its thermodynamic behaviour and stability by using MD simulation studies, employing "Academic version of Desmond" (Desmond 2020-3) Shaw Research, 2021). The System Builder tool was used for placing the apo protein or protein-ligand complex into an orthorhombic box, having a 10 Ð buffer region between the protein atoms and box sides, and filled with appropriate single point charge (SPC) water molecules. The system was neutralized with an adequate quantity of counter ions (Na þ and Cl -) and a fixed salt concentration of 0.15 M, representing the physiological monovalent ion concentration . Energy minimization was applied on the fully solvated and neutralized system to minimize any unfavourable interactions that may have occurred during the model building using the OPLS3e force field as the default settings in the Desmond panel. The temperature and pressure were set to 300 K and 1.01325 bar, respectively, for the isothermal-isobaric (NPT) ensemble (Kalibaeva et al., 2003;Martyna, 1994;Zrieq et al., 2021). For regulating the system temperature and pressure, the nose-Hoover thermostat and Martyna-Tobias-Klein algorithm was used. For short-range van der Waals and Coulomb interactions, a cut-off radius of 9.0 Ð was employed, while the Particle mesh Ewald method was used for the accurate and efficient evaluation of electrostatic interactions. A total of 100 nanosecond (ns) simulations were performed, with 1000 frames saved to the trajectory by recording the trajectory at 100 picosecond (ps) simulation time. Lastly, the simulation interaction diagram tool of the Desmond package was used to analyse the trajectory files . The outputs include the RMSD plots for the apo protein receptors and their respective complexes with the ligand, statistics of hydrogen bond, RMS fluctuation (RMSF), native ligand-protein contacts and radius of gyration (rGyr). The equations (i) and (ii) below were applied to compute the RMSD and RMSF (Shaw Research, 2021;Zrieq et al., 2021).

Dynamic Cross-correlation matrix (DCCM) and principal component analysis (PCA) of the MD simulation trajectories
In order to analyze the domain correlations, a dynamic crosscorrelation matrix (DCCM) was generated across all Ca-atoms for STOCK1N-45901 in complex with 7B3B, and 7BF6 during the 100 ns MD simulation. Protein-ligand confirmations and major global motions upon ligand binding was studied by principle component analysis (PCA). The DCCM and PCA were done using trj essential dynamics script of Schrodinger (Amadei et al., 1993;Bharadwaj et al., 2021).

Post simulation MM-GBSA binding free energy analysis
The thermal MM-GBSA.py script from the Prime module of the Schrodinger was used to carry out the post-simulation MM-GBSA analysis. Every 100th frame was collected from the last 50 ns of simulated trajectories from each MD trajectory, averaging over 6 frames, for lead compound binding free energy calculations. The overall ligand-protein complexes MM-GBSA binding free energy (DG Bind ) is calculated using the following equation (iii) (Zrieq et al., 2021).
where E ¼ prime energy.

Structure-activity relationship and prediction of antiviral-related activity
The SAR of the pharmacophoric sites of the finally selected natural product was studied about antiviral perspective using the available information from relevant literature. The contributions of the parent fragments, phenolic and aliphatic hydroxyl groups and the available positions for substitution to the interactive mechanisms of antiviral potentials were analysed relative to documented information (Annunziata et al., 2020;Bendary et al., 2013;Cramer et al., 2019;Harrold & Zavod, 2013;Lu et al., 2020;Stefaniu et al., 2020). Similarly, its bioactivities relevant to antiviral pharmacology were quantitatively predicted to further assess the possibility of proposed multi-target anti-SARS-CoV-2 mechanisms, using a server-based Prediction of Activity Spectra for Substances (PASS Online) (http://www.pharmaexpert.ru/passonline/index. php) by inputting the string of its SMILE format. The module applied a multilevel neighbourhood of atoms (MNA) structural descriptors of molecules to predict >4000 biological activities with %95% accuracy (leave-one-out cross-validation) (Filimonov et al., 2014), thus, reliable for predicting the activity of new compounds. The predicted results were in the numerical probability of active/inactive (pa/pi) respectively in SAR with the mechanisms of action and enzymatic pathways as antiviral, antiinfection, antiseptic, antitussive and peroxidase inhibitor.

Results and discussion
Active binding site cavities According to relevant previous studies Ni et al., 2021), the apo state crystallization of the SARS-CoV-2 nsp3 macrodomain constitutes two orthorhombic crystals of the enzyme, containing 3 and 5 molecules respectively contributed by Ile 131 and Phe 132, importantly for protein flexibility.
Other amino acid residues in the active binding pocket of the nsp3 macrodomain for pharmacological interactions include

Molecular docking
Molecular docking analyses were employed in series of SP and XP modules to hierarchically screen 102,174 protonated conformers produced from the ligand preparation step against the dual targets for SARS-CoV-2 replication-inclined polymerase activities, nsp3 and RdRp. The docking procedures were based on the searching algorithms of efficient heuristic ligand positioning and instant empirical scoring functions against the targets. The binding affinity is therefore basically estimated by the scoring function via the measurement of the noncovalent intermolecular interactions (enthalpic contributions) between the protein and the bound ligand after docking. Upon binding to small molecule ligand, the enzyme is generally believed to experience structural changes due to certain movement of the flexible side chain/backbone residues of the binding sites, allowing conformational complementarity of the enzyme with the shape and binding mode of the ligand. Thus the ligand selectively binds to a suitable conformational state of the enzyme active site (Du et al., 2016;Elekofehinti et al., 2021a). The protocols were validated by the redocking of the co-crystallized ligands onto the active sites, which produced poses with RMSD of <1.00 Ð against the protein structures. The resultant poses were comparable to the crystal structures retrieved from the RCSB PDB (Figure 4(A,B)), further validating the docking protocols as accurately reliable and precise (Castro-Alvarez et al., 2017;Ram ırez & Caballero, 2018). The XP glide scores of the finally screened 20 natural products with the most promising inhibitory interactions as well as the reference positive inhibitors, remdesivir and GS-441524 are presented (Table 1). From the results, it could be observed that the natural products STOCK1N-45901, -03804, -83408 and -08377 displayed consistently highest inhibitory interactions in terms of XP glide scores of -12.245, -11.701, -11.232 and -11.718 kcal/mol against the nsp3 and -8.864, -9.058, -8.574 and -9.116 kcal/mol against the nsp12, respectively. Interestingly, the 20 selected natural products all possess higher scores than remdesivir and GS-441524 references with XP values of -5.555 and -9.710 kcal/mol against the nsp3 and -6.168 and -4.545 kcal/mol against the nsp12, respectively. All the selected ligands and the references show a good glide energy model (emodel), supporting their strong inhibitory interactions.
The ligands interacted with the amino acid residues essential for the mitigation of the viral polymerase activities for replication, consistently with the references, remdesivir and GS-441524 as previously described. From the crystal structure of the nsp3 (PDB 7BF6), the binding interaction of GS-441524 as a small molecule inhibitor is similar to that of the ADP-ribose, a peptide target for removal during infection. Its NH group is sandwiched between amino acid residues Val 49 and Val 155, forming a H-bonding to Asp 22 (Figure 4(A)) (Ni et al., 2021). From the interactive binding pose obtained for the redocking of GS-441524 with the nsp3 macrodomain (Figure 4(B)), similar H-bonding interactions were observed through OH to Leu 126 and weak interactions to the backbone residues, Phe 156 and Asp 157. The four selected top-ranked ligands exhibited similar H-bonding interactions with the same set of amino acid residues and many other pharmacologically active ones through van der Waals non-bonding interactions, conferring high inhibitory interactions consistently with the XP docking scores. The details of the interactions are illustrated in 2D and 3D poses of the ligands at the active site of the enzyme ( Figure 5). The plasticity and flexibility of the protein active pocket reportedly constituted by some amino acid residues including Ile 131 and Phe 132 allows the binding of diverse inhibitors (Du et al., 2016;Marsh, 2015). These favour the interactive binding of the selected IBScreen natural products, STOCK1N-45901, -03804, -83408 and -08377 as promising inhibitors.
Similarly, from the previous binding interaction illustration (Figure 3(C)), remdesivir, a renowned inhibitor and co-crystallized ligand of the nsp12 structure interacted with amino acid residues, Asn 497, Lys 500 and Asp 845 as well as the RNA T-chain molecule T 8 all through H-bonding. It also formed p-p interactions with the RNA molecule T 9 and T 10 through its aromatic group, while other residues at the active site engaged in non-bonding van der Waals interactions for its inhibitory mechanisms of RdRp stalling. As shown in the generated 2D and 3D binding poses (Figure 6), among the four top-ranked natural product ligands with the highest XP scores, STOCK1N-45901 showed the more closely efficient ligand placement algorithms for similar biological action relative to remdesivir in terms of bonding interactions. It displayed strong interactions with amino acid residues Asn 496, Asn 497 and Lys 500, and RNA T 9, T 10 and T13 all through H-bonding, an indication for stronger RdRp stalling. Other selected natural products also displayed higher potentials than remdesivir in terms of binding interactions with the amino acid residues and RNA molecules. In addition to the H-bonding interactions commonly exhibited by the selected ligands as well as the references, STOCK1N-83408 displayed aromatic p-cation interaction with amino acid Lys 593 and aromatic p-p bonding with RNA molecules T 12 and T 12 along the chain A of the nsp12 active site. Similarly, STOCK1N-08377 formed a salt bridge with Lys 545 while the two reference inhibitors, remdesivir and GS-441524 also demonstrated aromatic p-p interactions with RNA molecules T 9, T 10, T 12 and T 13. Additionally, all the selected ligands including the reference inhibitors interacted with other pharmacologically essential amino acid residues including a side chain residue, Ser 861, an important target for rendering the nsp12 insensitive (stalled for RNA synthesis) upon binding to inhibitors (Kokic et al., 2021). The four natural product ligands, STOCK1N-45901, -03804, -83408 and -08377 consistently showed high XP scores and inhibitory binding poses against the dual targets relative to the co-crystallized known inhibitors, remdesivir and GS-441524. They are therefore adopted for further studies.

Prediction of pharmacokinetics, druggability and toxicological properties
The prediction of the physicochemical, pharmacokinetics, pharmacodynamics and toxicological properties are crucial for the identification of potential candidates with good delivery to targets for physiological expressions and excretion from the system without acute toxicity. The server-based pk-CSM module employed has been validated in alignment with several reported experimental results. Thus, it has ideally good accuracy and precision for predicting new chemical molecules under a drugability probe. The procedure is cost-saving and fast for identifying the most promising candidates subjectable to the time-consuming and expensive wet-lab modules for experimental validation (Pires et al., 2015). The predicted physicochemical properties of the natural product ligands selected from the docking studies and the reference inhibitors, GS-441524 and remdesivir are presented in Table  2. Only STOCK1N-45901 has a MW of <500 g/mol, others are heavier in weight. The Log P values of the four ligands were outside -0.7 to þ5.0, indicating poor lipophilicity. The ligands are highly flexible with numbers of rotatable bonds between 6-13. Again, only STOCK1N-45901 possesses <10 HBA, %5 HBD and TPSA 20-130 Å 2 . Advantageously, the MW of <500 g/mol and the number of rotatable bonds of 7 in the compound could influence >65% of its fractions for oral bioavailability of 20% and above in a whole rat system according to Veber's hypothesis (Veber et al., 2002). The displayed physicochemical parameters support that STOCK1N-45901 selectively has ideal flexibility, size, solubility, saturation and polarity, indicating ideal pharmacokinetics of lead-like prodrug candidate (Figure 7(A)) (Daina et al., 2017;Pires et al., 2015). Although, it is weakly lipophilic, however, this profile could be enhanced through a plausible structural modification using the SAR discussed later. The selected ligand satisfies the BDDC rule of 5 and druggability with insignificant violations and PAINS alert (Benet et al., 2016), as such, remains proficient for further assessment.
The pharmacokinetics and pharmacodynamics properties of the selected natural products, STOCK1N-45901, 03804, 83408 and 08377, and the references, GS441525 and remdesivir are shown in Table 2. Consistently with the physicochemical parameters, STOCK1N-45901 is distinctively predicted with a good solubility score and low permeability across the human intestinal epithelial Caco-2 cell. Although it shows lower intestinal absorption than STOCK1N-83408 and the standard drugs, however, it has a higher absorption potential than STOCK1N-03804 and STOCK1N-08377. Moreover, its modifiable structure and inertness for substitution could aid the enhancement for absorption in SAR. The four selected natural products mostly displayed no substrate expression for P-glycoprotein, a biological barrier by xenobiotics and extruding toxins, consistently with the remdesivir, except STOCK1N-83408. Again, STOCK1N-45901 demonstrated the highest fraction unbound in human for oral bioavailability compared to other ligands including the standard inhibitors. Several natural products generally possess low permeability across the BBB and the central nervous system (CNS), commonly to antiviral drugs due to high MW and low hydrophobicity (Bertrand et al., 2021), yet STOCK1N-45901 distinguishingly displayed high potential for BBB and CNS permeability. Interestingly, all the four selected natural products showed no substrate expression for metabolism-inclined cytochrome P450 isoenzymes, similarly to GS-441524, except STOCK1N-83408 for the metabolic CYP3A4. The ligand, STOCK1N-45901 consistently signals a good propensity for excretion from biological systems and low expressions across several pathways for toxicity including AMES, hERG I & II inhibition, hepatotoxicity, skin sensitisation, Tetrahymena Pyriformis, and minnow toxicity. Cumulatively, the predicted physicochemical, pharmacokinetics, pharmacodynamics and toxicological parameters for STOCK1N-45901 fall within the standards (Figure 7(A)) and distinctively support its promising potentials as a druggable candidate (Benet et al., 2016;Daina et al., 2017;Pires et al., 2015;Veber et al., 2002). The mutagenicity profile of the ligand read-across three similar molecules in SAR implies non-mutagenic (Figure 7(B)) (Gini et al., 2014;NaAllah et al., 2021), thus, it is ideally earmarked for further study.

Molecular dynamics
The conformational and thermodynamics stability of each target and in complex with the selected ligand was studied using MD simulations. This was employed to estimate the likely systemic biochemical interactive behaviour of the ligand virtually under specific conditions of physiological environments. Therefore, the promising lead natural product ligand, STOCK1N-45901 in complex with the SARS-CoV-2 nsp12 and nsp3 macrodomain was subjected to MD simulations to investigate its stability in the protein-ligand complex as well as main intermolecular interactions over the simulated trajectory. Additionally, the apo forms of both target structures were also simulated for comparative analysis. The RMSD and RMSF plots as well as the protein-ligand contacts were assessed from the resulting trajectories of the simulated complexes. The analysis of the RMSD plot (Figure 8(A)) indicates that the simulated system has equilibrated very well at the early period around 3 ns, then the fluctuations in the Ca atoms were consistently below 2.5 Ð during the entire simulated path of all complexes. The plot shows that after the initial period of RMSD fluctuations of proteins Ca atom due to  The natural product database annotations are abbreviated by removing the prefix, STOCK1N-.
equilibration, later the fluctuations in the Ca atoms were consistently throughout the simulated path of all complexes, an indication of good stability. Although, a slight fluctuation is not unexpected during the initial period but acquires stability throughout the rest of the simulation route. A system showing an average fluctuation of 1-3 Ð is usually Figure 5. 2D and 3D binding poses illustrating bonding and non-bonding interactions of selected natural products with amino acid residues at the active site of the nsp3 macrodomain (PDB 7BF6). Ball and stick models: Ligand (green), amino acid residue (red), RNA molecule (black). Interaction between ligand and residues are shown as hydrogen bonding (magenta arrow), p, cation (blue line), salt bridge (red-blue line) and p-p stacking (green line). The selected ligands showed more bonding interactions with active site amino acid residues and RNA molecules mostly through H-bonds and p-p interactions than remdesivir and GS-441524, indicating stronger binding affinities.
considered stable and deemed to be properly equilibrated in the case of globular proteins whereas elevated RMSD values indicate large conformational changes in protein structure over the progression of simulation Jamroz et al., 2013). Upon the MD simulation trajectories of 100 ns for the complexes of ligand, STOCK1N-45901 with the enzymes, nsp12 and nsp3 macrodomains (Figure 8(B)), slight fluctuations were observed within RMSD of 1-3 Å throughout, indicating stable systems (Jamroz et al., 2013;Patel et al., 2021). The rGyr of a ligand, comparable to its primary moment of inertia, is used to estimate its extended stability. Over a 100 ns simulation time, the rGyr plot (Figure 8(C)) also Figure 6. 2D and 3D binding poses illustrating bonding and non-bonding interactions of selected natural products with amino acid residues and RNA molecules at the active site of the nsp12 (PDB 7B3B). Ball and stick models: Ligand (green), amino acid residue (red), RNA molecule (black). Interaction between ligand and residues are shown as hydrogen bonding (magenta arrow), p, cation (blue line), salt bridge (red-blue line) and p-p stacking (green line). The selected ligands showed more bonding interactions with active site amino acid residues and RNA molecules mostly through H-bonds and p-p interactions than remdesivir and GS-441524, indicating stronger binding affinities.
shows that the lead compound, STOCK1N-45901 attained a good stability in each of the SARS CoV-2 nsp12 and nsp3 macrodomain binding pockets, with an average rGyr value of 3.99 ± 0.14 and 4.05 ± 0.19 Å, respectively (Table 3). More results of the 100 ns MD simulations performed for other lead compounds, STOCK1N-03804, STOCK1N-83408 and STOCK1N-08377 are available in supplementary information (supporting material Table S1and Figure S1). Investigation of hydrogen bonding pattern is used to evaluate the stability of protein structure and specificity of interaction which is associated with molecular recognition. Interactions of hydrogen bonds between ligand (STOCK1N-45901) and dual target proteins of the SARS-CoV-2 through the 100 ns of MD simulations were evaluated (Figure 9(A,B)). An average of intermolecular hydrogen bonds of STOCK1N-45901 in complex with nsp12 (PDB 7B3B) and nsp3 macrodomain (PDB 7BF6) were found to be 772 and 153, respectively, which confirmed the strong inhibition of targeted protein.
The RMSF plots (Figure 10(A)) are used to examines the regions of protein structure that are fluctuating from their mean structure. The a-helical, b-strand and loop regions are depicted in red, blue and white backgrounds respectively. In contrast to the loop area, a-helical and b-strands regions of the proteins are stiffer than the unstructured regions. It could be inferred that the lead compound, STOCK1N-45901 was tightly bound within the respective cavity of the SARS-CoV-2 nsp12 and nsp3 macrodomain protein binding pocket, because the active site and main chain atoms moderately fluctuated mostly within 1-3 Ð, indicating a good stability . Observably, the ligand contacted 20 amino acids of nsp12 protein,namely,Ile494,Val495,Asn496,Asn497,Leu498,Lys500,Lys545,Tyr546,Arg569,Gln573,Lys577,Ile589,Gly590,Thr591,Ser592,Trp598,Ala688,Tyr689,Leu758 and Cys813. In nsp3 macrodomain protein, it interacted with 28 amino acids including Ala21, Ile23, Asn37, Ala38, Ala39, Asn40,Val41,Lys44,His45,Gly46,Gly47,Gly48,Val49,Ala50,Ala52,Gly97,Pro98,Asn99,Pro125,Leu126,Ser128,Ala129,Ile131,Phe132,Ala154,Phe156,Asp157,and Leu160. Each interaction of residue is represented by a vertical bar in green color. All these interacted residues have a RMSF value of less than 1.2 Ð, which indicates that reported compound, STOCK1N-45901 is tightly bound to the dual target proteins of the SARS-CoV-2. The RMSF plot revealed an insignificant fluctuation, consistently with the RMSD and rGyr plots. The Statistics of hydrogen bonding chart ( Figure  10(B)) also indicates the bonding interactions made by the ligand with amino acid residues of the proteins during simulations, contributing to the high binding affinity of the ligand to the enzyme targets. The illustration of native contact made with the active amino acid residues (Figure 10(C)) showing the interactions of the pharmacophoric sites of the ligand that result in strong inhibitory binding to the receptors during simulations. Figure 10(B) also depicts the protein-ligand complex binding interaction histogram during the simulation time. In nsp12-STOCK1N-45901 complex, the amino acids Ile494, Val495, Asn496, Asn497 and Lys500 contributed primarily to the both hydrogen bonding and water mediated hydrogen bond interactions (5%-65%), while Ile589 Ala688 and Leu758 contributed mainly to the hydrophobic interactions (5%-19%). Protein ligand contact analysis of STOCK1N-45901 in complex with nsp12 shows that the residues Asn496 and Asn497 exhibited 21% and 11% hydrogen bond interactions, respectively, with the terminal hydroxyl group of lead compound (Figure 10(C)). For STOCK1N-45901-nsp3 macrodomain complex hydrogen bond interaction is seen with Asn40, Lys44, Gly46, Gly47, Leu126, Ser128, Ala129, and Ala154. Water bridges with these residues were also found as key interaction bonds. Also, Major hydrophobic interactions were seen with Ala38, Val49, Phe132 and Phe156. nsp3 macrodomain amino acid residue ser128 exhibited bidentate hydrogen bond interaction with STOCK1N-45901, similarly Lys44 and Leu126 also interacted with STOCK1N-45901through one hydrogen bonding.

Post-simulation MMGBSA analysis
The post simulation MM-GBSA analysis of binding free energy (DG Bind ) calculation was carried out with the generation of 1-1000 frames having a 1-step sampling size. A total of 1000 frames were processed and analysed throughout the 100 ns MD data of the lead compound, STOCK1N-45901 in complex with the nsp12 and nsp3 macrodomain revealed by the dynamics studies ( Figure 11). The computed per frame post simulation-MMGBSA binding free energy for the protein ligand complexes is depicted in supporting material Tables S2 and S3. The calculated average DG Bind of the product, STOCK1N-45901 in complex with the SARS CoV-2 nsp12 and nsp3 macrodomain was found -60.70 ± 7.73 kcal/mol, and -55.15 ± 4.75 kcal/mol, respectively, indicating strong bindings (Table 3) (Halgren et al., 2004).

Dynamic cross-correlation matrix (DCCM) and principal component analysis (PCA) of the MD simulation trajectories
The DCCM was developed to investigate the correlated mobility of structural domains in order to achieve a stable conformation of nsp12 and nsp3 macrodomain protein after STOCK1N-45901 binding from MD trajectories. This is calculated through the calculation of all pairwise correlation coefficients. The correlation values are calculated between -1 and þ1, where colors ranging from red to blue reflect the    intensity of associated motion between residues, with blue colors indicating positive (þ1/complete) correlation, white indicating no correlation, and red indicating negative (-1/ complete anti) correlation. The correlation scores on the central mean line (blue) revealed four and three separate blocks with a high correlation of amino acid movement in nsp12 and nsp3 macrodomain protein, respectively. In nsp12, the domain D1 have highest cross-correlation between residues comprised of 36-152 residues conforming into two distinct b-sheets, one helix and an extended loop (red) whereas, D2 domain having the 180-223 residues conforming into more flexible extended loops and one b-sheets (blue). The D3 and D4 domains conform to 380-460 (magenta) and 550-680 (yellow) residues, respectively ( Figure 12(A)). When comparing the DCCM of nsp3 macrodomain protein, it was observed that the residues in the range of 0-20 in D1(red), 36-84 in D2(blue), and  residues in the range of 105-145 in D3(pink) show highly correlated directions of dynamical movement (Figure 12(B)). The DCCM analysis confirms good agreement with the RMSF of the C-a backbone of nsp12 and nsp3 macrodomain proteins in the STOCK1N-45901 bound state, with moderate to fewer fluctuations of respective amino acid residues. On the trajectories generated by MD simulations, we also did PCA to investigate conformational distribution during the simulation time and investigate large-scale collective motions of the protein in protein-ligand complexes. For all of these complexes, projection of the protein motion in phase space along the PC1 and PC2 revealed a uniform distribution of conformations throughout the simulations (Figure 13).

SAR and prediction of antiviral-related activity
The SAR study of the selected natural product ligand, STOCK1N-45901 shows the relevance of its pharmacophoric fragments to antiviral activity using the available information from relevant literature (Annunziata et al., 2020;Bendary et al., 2013;Cramer et al., 2019;Lu et al., 2020;Stefaniu et al., 2020). The entire structure ( Figure 14) allows a compact molecular spatial fit into active sites of biological targets. The parent, 2-hydroxy benzoate fragment is renowned for various pharmacological activities including anticancer, antifeedant, antifungal, antiinfection, anti-inflammation, antioxidant, L-ascorbate inhibitor and promoter of exfoliation. The phenolic and aliphatic hydroxyl groups are prominent flexible positions for substitution for enhancing the ligandbased mechanisms of antiviral potentials such as lipophilicity and binding affinity. The ortho (position 2/6) phenolic OH is more favoured for activity-related intermolecular H-bonding than meta (position 3/5) and para (position 4) while the lower bond dissociation energy supports antioxidant effects. Additionally, the presence of polyols commonly in natural products favour pharmacodynamic effects. Although, many aliphatic hydroxyl groups in original structures depict a propensity for high desolvation penalty which could affect binding affinity, commonly to naturally derived prodrugs. Moreover, metabolic pathways affect the functions of OH groups depending on their stereochemistry. For instance, primary and phenolic OH groups have a strong tendency for rapid metabolism and production of toxic metabolites. However, the flexible chemistry could aid beneficial modulation of the shortfalls with appropriate substitution by more suitable fragments (Annunziata et al., 2020;Bendary et al., 2013;Cramer et al., 2019;Harrold & Zavod, 2013;Lu et al., 2020;Stefaniu et al., 2020).
Similarly, the bioactivities of the ligand, STOCK1N-45901 relevant to antiviral pharmacology were quantitatively predicted to further assess the possibility of proposed multi-target anti-SARS-CoV-2 mechanisms, using an MNA molecular descriptor incorporated in PASS Online. The predicted results in SAR with the relevant mechanisms of action and enzymatic pathways are 0.656, 0.713, 0.635, 0.827, 0.655 and 0.815 for the probability of activity and 0.004, 0.007, 0.025, 0.004, 0.005 and 0.004 for the probability of inactivity as antiviral, antiinfection, anti-inflammation, antiseptic, antitussive and peroxidase inhibitor respectively. The accuracy of the module has been validated %95% with the chance for its attainment when the activity predicted for a molecule is !5, as obtained in this study. Previous studies have used the module to predict biological activities of some natural products, the results of which were aligned with experimental outcomes (Filimonov et al., 2014;Grienke et al., 2015). These support the reliability of the predictions and pose the finally retrieved ligand, STOCK1N-45901 as a promising active multi-target directed antiviral candidate for further transformational assessment against SARS-CoV-2.

Conclusion
The nsp12 (RdRp) and nsp3 macrodomain enzymes are strongly implicated in the virulence regulation of host immune and polymerase activities for the replication of SARS-CoV-2. These make them plausible therapeutic targets for effective mitigation of SARS-CoV-2 infectivity. However, remdesivir remains the only successful FDA-approved small molecule mechanistic inhibitor of the RdRp in clinical conditions of the COVID-19, while none has been reported for the nsp3 macrodomain. The preliminary application of molecular modelling strategies involving SP and XP ligand-receptor docking screening in this study has afforded the identification of four ligands out of the 69,067 natural products retrieved from the IBScreen database with the highest inhibitory potency against the dual targets. Further assessments through the predictions of physicochemical, pharmacokinetics, pharmacodynamics, druggability, toxicology and mutagenicity parameters distinguished the ligand, STOCK1N-45901 as relatively outstanding. The ligand, upon subjection to 100 ns MD simulation trajectories displayed interesting conformational and thermodynamic stabilities in complex with the protein targets. In a SAR, its molecular structure favours biological activities relevant to the proposed anti-SARS-CoV-2 pharmacology, in alignment with documented literature and PASS Online activity predictions. Interestingly, the molecular structure supports flexible chemistry for modulating the desired biological activities and suitability for effective one-drug-multiple-target anti-SARS-CoV-2 pharmacology. Although, rigorous wet-lab biochemical analyses are  required for more adequate confirmation of the screened candidate as a bona fide inhibitor of the dual targets of SARS-CoV-2. Moreover, the final compound recommended for further study contains phenol that may undergo tautomerisation in its physical state through the transfer of phenolic H to the C at the ortho position, commonly to %25% of database compounds. Such investigation is important for the experimental validation of the compounds as it may affect the physical and biochemical properties including hydrophobicity and binding affinity. As such, further experimental investigations on the finally retrieved compound deserves tautomerisation and identification of the specific tautomer for the proposed activity. However, the results in this study build confidence for embarking on such expensive, time-demanding and rigorous evaluations and pose the natural product, STOCK1N-45901 (2,3,4,5,6-pentahydroxyhexyl o-hydroxybenzoate) from the IBScreen screening database as a promising inhibitor of the dual RdRp and nsp3 of SARS-CoV-2 amenable for further study.