Discovery of Bispecific Lead Compounds from Azadirachta indica against ZIKA NS2B-NS3 Protease and NS5 RNA Dependent RNA Polymerase Using Molecular Simulations

Zika virus (ZIKV) has been characterized as one of many potential pathogens and placed under future epidemic outbreaks by the WHO. However, a lack of potential therapeutics can result in an uncontrolled pandemic as with other human pandemic viruses. Therefore, prioritized effective therapeutics development has been recommended against ZIKV. In this context, the present study adopted a strategy to explore the lead compounds from Azadirachta indica against ZIKV via concurrent inhibition of the NS2B-NS3 protease (ZIKVpro) and NS5 RNA dependent RNA polymerase (ZIKVRdRp) proteins using molecular simulations. Initially, structure-based virtual screening of 44 bioflavonoids reported in Azadirachta indica against the crystal structures of targeted ZIKV proteins resulted in the identification of the top four common bioflavonoids, viz. Rutin, Nicotiflorin, Isoquercitrin, and Hyperoside. These compounds showed substantial docking energy (−7.9 to −11.01 kcal/mol) and intermolecular interactions with essential residues of ZIKVpro (B:His51, B:Asp75, and B:Ser135) and ZIKVRdRp (Asp540, Ile799, and Asp665) by comparison to the reference compounds, O7N inhibitor (ZIKVpro) and Sofosbuvir inhibitor (ZIKVRdRp). Besides, long interval molecular dynamics simulation (500 ns) on the selected docked poses reveals stability of the respective docked poses contributed by intermolecular hydrogen bonds and hydrophobic interactions. The predicted complex stability was further supported by calculated end-point binding free energy using molecular mechanics generalized born surface area (MM/GBSA) method. Consequently, the identified common bioflavonoids are recommended as promising therapeutic inhibitors of ZIKVpro and ZIKVRdRp against ZIKV for further experimental assessment.


Introduction
Zika virus (ZIKV) was first isolated in 1947 from Zika forest, Uganda, East Africa [1], and remained unnoticed for almost 60 years. In 2007, this virus caught everyone's attention following the first-ever ZIKV epidemic outbreak in Yap Island, Federated States of Micronesia, where 59 predictable and 49 confirmed ZIKV cases were reported [2]. Since then, ZIKV has caused several epidemics outside African countries in the last ten years, including the 2013-2014 outbreak in French Polynesia, infecting around 28,000 people [3,4]. Subsequently, in 2015, suspected outbreak of ZIKV in Brazil was estimated to infect 440,000 to 1,300,000 individuals [5], while microcephaly and other neurological disorders were also observed in approximately 7000 infected individuals [6,7]. In 2016, several cases of ZIKV infection were observed in females from the United States of America (USA) who had never travelled to the countries affected with this virus, but their male partners did [8]; the presence of ZIKV in their semen confirmed that it could also be transmitted through sexual contact [9]. Notably, similar to all of the flaviviruses, ZIKV is also primarily transmitted through Aedes aegypti mosquitoes [10]. However, the transfer of ZIKA infection through sexual transmission [8] and the vertical transmission from mother to the fetus [11][12][13] marks this virus as a global health concern. Currently, no therapeutics or treatments are available for the ZIKV infection. As a consequence, ZIKA is posing a serious threat to humans globally. Thus, this raises a demand for the development of potential therapeutics to control the ZIKA epidemic and associated neurological disorders.
ZIKV is a vector-borne envelop flavivirus and encloses a 10.8 kb positive sense, singlestranded, RNA (+ssRNA) genome, which contains a single open reading frame (ORF) for the translation of a single polyprotein of 3419 amino acids [14]. Genome replication plays a central role in the viral pathogenesis. Thus, after infection, ZIKV polyprotein is processed into three structural proteins (SPs): pre-membrane (prM), envelope (E), and capsid (C) proteins, and seven nonstructural proteins (NSPs): NS1, NS2A, NS2B, NS3, NS4A, NS4B, and NS5, via the proteolytic activity of both ZIKA and host proteinases ( Figure 1) [15]. The structural proteins provide the protection to the newly synthesized viral genome by forming an inner layer of capsid proteins while the precursor membrane (prM) protein and an envelope (E) protein contribute to the formation of the virion surface. During maturation, the prM protein is proteolytically cleaved into pr-subunit and M-subunit by the catalytic activity of host furin protease in the trans-Golgi network (TGN). This event results in the formation and release of fully matured ZIKV with E and M protein on its outer envelope from the host cell [16]. Genome replication is the crux of viral pathogenesis, and in the case of ZIKV, the NSPs interact to form a replication complex that provides a site for the synthesis of RNA genome of the viral genomic RNA. Among all of the NSPs, NS2B-NS3 protease (ZIKV pro ) and NS5 RNA dependent RNA polymerase (ZIKV RdRp ) are the vital factors in ZIKV pathogenesis, as the former one is involved in the hydrolysis and maturation of the flavivirus polyprotein whereas the latter one has polymerase activity, which is necessary for the viral replication process [17].
The ZIKV serine protease (ZIKV pro ), a heterodimeric complex, is consists of a membrane protein NS2B bound with~70 kDa NS3 protein at the N-terminal region [18][19][20]. The NS3 protein has protease and helicase domains at the N-terminal and the C-terminal, respectively. However, despite lacking any enzymatic activity, NS2B plays a crucial role in the folding of NS3 protein [21][22][23]-acted as a co-factor for the activity of NS3 protein [24] and holds it near the cell membrane, which is essential for its proteolytic activity and viral replication [25][26][27][28][29]. Thus, in NS2B-NS3 protease (ZIKV pro ), the substrate binding and catalyzing active site of the NS3 protease domain is enfolded by the NS2B protein. Herein, a stretch of forty amino acid residues located at the C-terminal region of NS2B interact with the N-terminal protease domain of NS3 protein [30][31][32][33][34][35], which results in the formation of a catalytic triad (His 51 , Asp 75 , and Ser 135 residues). These catalytic residues are required for the proteolytic activity by the virus to release the functional NSPs in the cytosolic side of the host endoplasmic reticulum (ER), which further participate in viral replication ( Figure 1) [36,37]. Due to this crucial role of ZIKV pro in the life cycle ZIKV and the lack of any a homolog in human cells, ZIKV pro is considered as a promising target for anti-ZIKV drug development.
Moreover, NS5 is the largest and highly conserved protein among flaviviruses, including ZIKV, which showed essential function in the viral genome replication via the Cterminal NS5 RdRp (ZIKV RdRp ) domain. Therefore, targeting the ZIKV RdRp domain has been considered as a precise therapeutic strategy against ZIKV [38][39][40]. The structural analysis revealed the right-hand-shaped conformation for the ZIKV RdRp domain, holding three main domains: fingers, palm, and thumb, where the finger and thumb domains intersect to form an active region with a central catalytic pocket formed by the palm domain. Herein, the amino acids ranging from 321 to 488 and 542 to 608 comprise the finger domain, 484 to 541 and 609 to 714 comprise the palm domain, and 715 to 903 comprise the These catalytic residues are required for the proteolytic activity by the virus to release the functional NSPs in the cytosolic side of the host endoplasmic reticulum (ER), which further participate in viral replication ( Figure 1) [36,37]. Due to this crucial role of ZIKV pro in the life cycle ZIKV and the lack of any a homolog in human cells, ZIKV pro is considered as a promising target for anti-ZIKV drug development.
Moreover, NS5 is the largest and highly conserved protein among flaviviruses, including ZIKV, which showed essential function in the viral genome replication via the C-terminal NS5 RdRp (ZIKV RdRp ) domain. Therefore, targeting the ZIKV RdRp domain has been considered as a precise therapeutic strategy against ZIKV [38][39][40]. The structural analysis revealed the right-hand-shaped conformation for the ZIKV RdRp domain, holding three main domains: fingers, palm, and thumb, where the finger and thumb domains intersect to form an active region with a central catalytic pocket formed by the palm domain. Herein, the amino acids ranging from 321 to 488 and 542 to 608 comprise the finger domain,  484 to 541 and 609 to 714 comprise the palm domain, and 715 to 903 comprise the thumb  domain, where Asp 540 (palm domain) and Ile 799 and Asp 665 residues (thumb domain) were identified to form the catalytic site of ZIKV RdRp and crucial for the interaction with the ligands (Figure 1) [41,42]. In addition, NS5 protein also carries methyltransferase (MTase) activity at the N-terminal end, which is required for the 5 capping of newly synthesized viral mRNA [43].
In 2016, after the WHO announced the ZIKV outbreak with a global health emergency, various therapeutics, such as orthosteric inhibitors [44,45], allosteric inhibitors [46][47][48], ZIKV pro inhibitors [49,50], ZIKV RdRp inhibitors [51][52][53][54], and a few inhibitors with unknown molecular targets [55][56][57], were reported as a treatment for the ZIKV infection. However, only one compound, viz. novobiocin, was noted for considerable in vivo inhibitory effect against the ZIKV infection [58]. Therefore, in the war against ZIKV, a comprehensive blueprint is needed to develop a promising anti-ZIKV therapeutics. In this context, among the various therapeutic designing methods, a multitargeted approach has been suggested as an aspiring strategy where the most appealing targets for the ZIKV are ZIKV pro and ZIKV RdRp domains.
In the last two decades, the development of multitargeted drugs has been taken into preference due to their major advantages, due to the lower risk for drug interactions and improved drug compliance in patients [59][60][61]. In this context, the present study opted the multitargeting approach against ZIKV by identification of potent bioflavonoids as ZIKV pro and ZIKV RdRp domain putative inhibitors from the Azadirachta indica, popularly known as Neem, which is well established to possess antibacterial, antifungal, and antiviral activity [62]. Of note, the various parts of the A. indica, such as leaves, flowers, bark, seeds, and roots, are used in several therapies and treatments for the infectious and noninfectious diseases in Asian and African countries since time immemorial. In recent years, medicinal plants, including A. indica, are the foremost choice in finding a cure against numerous diseases due to their least toxicity, unique chemistry of secondary metabolites, and a long-term resource with constant mass production [63,64]. Therefore, in this study, 44 bioflavonoids reported from Azadirachta indica, were computationally screened in the active pocket of ZIKV pro and ZIKV RdRp domains to identify bispecific potent inhibitors with substantial binding affinity and stability for the drug development against ZIKV infection.

Receptors and Bioflavonoids
The three-dimensional (3D) crystal structure of ZIKV NS2B-NS3 protease (ZIKV pro , PDB ID: 6Y3B [65]) and ZIKV NS5 RNA-dependent RNA polymerase domain (ZIKV RdRp , PDB ID: 6LD2 [52]) of 1.59 and 1.40 Å resolutions, respectively were downloaded from the protein data bank (PDB) database (https://www.rcsb.org/) [66]. The selected proteins as receptors were preprocessed by assigning bond order and addition of polar hydrogen atoms using the default parameters in Protein preparation wizard of the Maestro-Schrödinger suite [67]. Herein, protein structures were treated for protonation of residues using the PROPKA program at pH 7.0, followed by the restrained minimization using Optimized Potentials for Liquid Simulations 3e (OPLS3e) force field under default parameters.
To identify the bioflavonoids from Azadirachta indica (Neem plant) as putative inhibitors of ZIKV pro and ZIKV RdRp , a small library of known 44 bioflavonoids was prepared by exploring the documented research articles (Table S1). The three-dimensional conformers of all bioflavonoids were retrieved from the PubChem database (https://pubchem.ncbi. nlm.nih.gov/) [68] and treated as ligand for the computational analysis against targeted ZIKV proteins. Briefly, the ligand library was prepared under default parameters using LigPrep panel in Schrödinger suite (Schrödinger Release 2018-3: LigPrep, Schrödinger, LLC, New York, NY, USA, 2018), where ligand tautomeric conformations were generated using EPIK state penalty at pH 7.0 ± 2.0 with OPLS3e force field.

Structure-Based Virtual Screening and ADMET Analysis
In the initial stage of drug discovery, structure-based virtual screening (SBVS) plays a crucial role in identifying the novel bioactive molecules as potent ligand against the threedimensional structure of a certain biological targets obtained through X-ray diffraction, NMR, Cryo-electron microscopy, or homology model. SBVS is a computational technique which attempts to predict the putative conformations between the receptor and ligand for complex formation and uses the non-covalent interactions-based scoring function to mark the stability of calculated receptor-ligand complexes.

Redocking and Intermolecular Interaction Profiling
The top common compounds collected from SBVS against ZIKV proteins, i.e., ZIKV pro and ZIKV RdRp , and respective reference compounds, i.e., O7N inhibitor (native co-crystalized ligand) for ZIKV pro and Sofosbuvir inhibitor for ZIKV RdRp (previously reported nucleoside inhibitor of ZIKV RdRp ) [72], were redocked in the selected respective binding pockets of viral proteins under default parameters using Glide XP module of Schrödinger suite (Schrödinger Release 2018-3: Glide, Schrödinger, LLC, New York, NY, USA, 2018). All of the docked poses were studied for the intermolecular interactions under the default parameters in the Maestro v12.9 package, and both 3D and 2D interaction diagrams were prepared using the free academic version of the Maestro v12.9 package (Schrödinger Release 2021-3: Maestro, Schrödinger, LLC, New York, NY, USA, 2021).

Molecular Dynamics Simulation Analysis
Dynamic stability and the intermolecular interactions profiling of the selected proteinligand complexes were analyzed through the molecular dynamics (MD) simulation using free academic Desmond-maestro 2020-4 package [73,74]. Initially, each docked complex was placed in a 10 Å × 10 Å × 10 Å orthorhombic box amended with explicit (TIP4P: transferable intermolecular potential 4 point) solvent using a system builder module. Following, the complete simulation system was amended with 0.15 M salt to mimic the physiological conditions and neutralized using counter sodium and chlorine ions while placed at 20 Å distance from the docked ligand in the binding pocket of the receptor. Later, the simulation system was minimized under default parameters using a minimization tool and subjected to 500 ns simulation under OPLS-2005 force field at 300 K temperature and 1.01325 bar pressure with default parameters using Molecular dynamics simulation tool of free academic Desmond-maestro 2020-4 [73,74]. At last, the MD simulation trajectory of each protein-ligand complex was analyzed for the stability and intermolecular interactions as a function of 500 ns interval by simulation interaction diagram (SID) module in the free academic Desmond-maestro 2020-4 suite [73,74].
The following Equation (1) calculated the root mean square deviation (RMSD) values for the protein alpha carbon (Cα) atoms and the ligand heavy atoms with respect to protein (Cα) in each frame during the 500 ns simulation trajectory to measure the average deviation that occurred in the protein and ligand for the respective docked complex in reference to their respective initial poses [75].
While calculating RMSD values, N represents the number of atoms selected; t ref is defined as reference time at zero interval; r i denotes the position of the atoms under evaluation in frame x followed by the superimposition on the reference frame r i at time interval t x .
Moreover, root mean square fluctuation (RMSF) values were also calculated for characterizing the local fluctuations at residue and atomic level in protein structure and ligand molecule, respectively. The following equation (2) expresses the local fluctuation in the simulation trajectory [75].
While calculating RMSF, T denotes the simulation interval for which the RMSF is calculated, t ref denotes the reference time, r i denotes the atom position in reference time t ref and r i denotes atom position at the time following superimposition on the reference frame.

Endpoint Free Binding Energy Calculation
Molecular mechanics/generalized Born Surface area (MM/GBSA) analysis was performed to calculate the mean binding free energy on the extracted poses from the last 10 ns interval (at 10 ps step) of respective MD simulation trajectory under default parameters with OPLS-2005 force field in the Prime MMGBSA module in Schrödinger suite (Schrödinger Release 2018-3: Prime, Schrödinger, LLC, New York, NY, USA, 2018). Herein, all of the solvent molecules and ions were deleted from the extracted poses, and the binding free energy (∆G) was calculated using the following Equation (3).
Where, ∆G Bind = Binding free energy, ∆G Complex (minimized) = Free energy of the complex, ∆G Receptor (minimized) = Free energy for the receptor, and ∆G Ligand (minimized) = Free energy for the ligand.

Structure-Based Virtual Screening
The primary goal of this study was to find the common compounds from a natural source that can inhibit both the ZIKV pro and ZIKV RdRp proteins for the treatment of ZIKV infection. Thus, a small library of 44 bioflavonoids (Table S1) belonging to the A. indica was used in SBVS against the selected binding pocket of ZIKV pro and ZIKV RdRp . This resultant in the collection of 21 compounds with docking scores between −2.0 to −11.01 kcal/mol against the selected viral proteins (Tables S2 and S3). Following, based on their docking scores, only the top four common bioflavonoids, i.e., Rutin, Nicotiflorin, Isoquercitrin, and Hyperoside, were marked as bispecific inhibitors for further redocking and intermolecular interaction (IMI) analysis by comparison to the reference compounds of ZIKV pro and ZIKV RdRp , i.e., O7N inhibitor and Sofosbuvir inhibitor, respectively ( Figure 2). Herein, the selected four bioactive bioflavonoids showed significant docking scores between −8 to −11.1 kcal/mol with the target proteins, i.e., ZIKV pro and ZIKV RdRp domain (Tables S2 and S3). Interestingly, all of the identified bioflavonoids were previously reported to have medicinal and therapeutic properties; for instance, Rutin and Isoquercitrin were documented for antiviral, anticancer, and antidiabetic activities [76][77][78][79][80][81][82][83], Nicotiflorin was described to inhibit SARS-CoV-2 M pro and dengue virus NS2B-NS3 protease [75,84,85], and Hyperoside was also reported to have anticancer activity [86,87].

Redocking and Intermolecular Interaction Analysis
Redocking is a mandatory analysis after SBVS calculation to assure that the compounds identified and selected through virtual screening have high affinity with the active site residues of the binding pocket since the algorithms of SBVS are fast and, therefore, their accuracy level is comparatively low than the docking scoring methods [88]. Thus, a stringent XP docking protocol was adopted in the redocking of the selected poses, and the most satisfactory binding poses with substantial binding scores and interactions with the essential residues in the viral proteins, i.e., ZIKV pro and ZIKV RdRp , were extracted for further analysis. Herein, the redocked complexes of ZIKV pro with Rutin, Nicotiflorin, Isoquercitrin, and Hyperoside were noted for docking energy of −10.61, −9.95, −8.63, and −8.37 kcal/mol, respectively (Table 1). Likewise, Rutin, Nicotiflorin, Isoquercitrin, and Hyperoside compounds docked with ZIKV RdRp showed higher docking scores of −11.01, −10.56, −8.84, and −7.87 kcal/mol, respectively (Table 1). Interestingly, all four bioactive bioflavonoids, i.e., Rutin, Nicotiflorin, Isoquercitrin, and Hyperoside, demonstrated higher redocking scores (−7.8 to 11.01 kcal/mol) with both the viral target proteins by comparison to the respective reference inhibitors, viz. O7N inhibitor for ZIKV pro (−6.629 kcal/mol) and Sofosbuvir inhibitor for ZIKV RdRp (−6.033 kcal/mol). Therefore, the redocking results

Redocking and Intermolecular Interaction Analysis
Redocking is a mandatory analysis after SBVS calculation to assure that the compounds identified and selected through virtual screening have high affinity with the active site residues of the binding pocket since the algorithms of SBVS are fast and, therefore, their accuracy level is comparatively low than the docking scoring methods [88]. Thus, a stringent XP docking protocol was adopted in the redocking of the selected poses, and the most satisfactory binding poses with substantial binding scores and interactions with the essential residues in the viral proteins, i.e., ZIKV pro and ZIKV RdRp , were extracted for further analysis. Herein, the redocked complexes of ZIKV pro with Rutin, Nicotiflorin, Isoquercitrin, and Hyperoside were noted for docking energy of −10.61, −9.95, −8.63, and −8.37 kcal/mol, respectively (Table 1). Likewise, Rutin, Nicotiflorin, Isoquercitrin, and Hyperoside compounds docked with ZIKV RdRp showed higher docking scores of −11.01, −10.56, −8.84, and −7.87 kcal/mol, respectively (Table 1). Interestingly, all four bioactive bioflavonoids, i.e., Rutin, Nicotiflorin, Isoquercitrin, and Hyperoside, demonstrated higher redocking scores (−7.8 to 11.01 kcal/mol) with both the viral target proteins by comparison to the respective reference inhibitors, viz. O7N inhibitor for ZIKV pro (−6.629 kcal/mol) and Sofosbuvir inhibitor for ZIKV RdRp (−6.033 kcal/mol). Therefore, the redocking results concluded that each of the selected conformations of the docked bioflavonoids have established a considerable binding affinity with the binding pocket of selected viral targets, i.e., ZIKV pro and ZIKV RdRp , and can considered for computational analysis. Table 1. Redocking score and intermolecular interactions noted for the screened compounds with the viral proteins, i.e., ZIKV pro and ZIKV RdRp , within 4 Å around the docked ligand in the respective binding pockets.
Intermolecular interaction (IMI) analysis is essential to understand the mode of molecular contact formation between the docked ligands and the target proteins. Herein, each docked bioflavonoid (Rutin, Nicotiflorin, Isoquercitrin, and Hyperoside) was observed for the formation of hydrogen bond (H-bond), hydrophobic, and π-π interactions with the essential residues of target proteins (ZIKV pro and ZIKV RdRp ) ( Table 1 and Table S4 and Figures 3 and 4). In details, the docked complex of ZIKV pro −Rutin was observed for the formation of four H-bonds via B:Val 36 , A:Ser 81 , B:Asn 152 , and B:Gly 153 residues, and two π-cation stacking interactions with B:His 51 residue (Figure 3). Also, ZIKV pro −Nicotiflorin  (Figure 3). Additionally, all ZIKV pro −bioflavonoids docked complexes were identified for intermolecular hydrophobic, polar, negative, positive, and glycine interactions (Table 1 and Table S4, Figure 3).    Furthermore, the docked complex of ZIKV RdRp −Rutin was also observed for the formation of seven H-bonds with Glu 419 , Gly 604 , Asp 666 , Ser 798 , and Ile 799 residues and two π-π stacking interactions with Trp 797 residue (Figure 4). Whereas the ZIKV RdRp −Nicotiflorin docked complex exhibits the formation of six H-bonds with Asp 540 , Trp 539 , Asp 665 , Asp 666 , and Cys 711 residues while ZIKV RdRp −Isoquercitrin complex showed only five H-bonds with Asp 540 , Asp 665 , Asp 666 , Cys 711 , and Ile 799 residues (Figure 4). Compared to other complexes, ZIKV RdRp -Hyperoside docked complex included only two amino acid residues (Asp 540 , and Asp 666 ) to form four hydrogen bonds (Figure 4). Additionally, hydrophobic, polar, positive, negative, and glycine interactions were also observed in all of the ZIKV RdRp −Ligand docked complexes (Figure 4, Table 1 and Table S4). Notably, a similar intermolecular interaction profile was observed in the reference docked complexes, i.e., ZIKV pro −O7N (Figure 3) and ZIKV RdRp −Sofosbuvir (Figure 4, Table 1 and Table S4). Collectively, the analysis of interaction profiles of all of the docked poses advocates the identified bispecific bioflavonoids for occupying similar active regions in targeted viral protein with higher binding energy by comparison to the respective reference inhibitors.

ADMET and Drug-Likeliness Analysis
In the field of drug discovery, the compounds or molecules proposed as a drug candidate must carry high biological activity and no or least toxicity. Therefore, a few critical pharmacological parameters, such as absorption, distribution, metabolism, excretion, and toxicity (ADMET parameters) along with the pharmacokinetics, has been suggested for the validation on every proposed drug candidate. Early assessments of such parameters in the initial phase of drug discovery are essential to understand and avoid drug molecules' pharmacokinetics-related failure during clinical trials [89]. Thus, to analyze the pharmacokinetic properties and drug-likeness, all of the screened bioactive bioflavonoids, i.e., Rutin, Nicotiflorin, Isoquercitrin, and Hyperoside, as well as the reference compounds, i.e., O7N inhibitor (for ZIKV pro ) and Sofosbuvir inhibitor (for ZIKV RdRp ) (Figure 2), were uploaded on the SwissADME and admetSAR online servers for the assessment of AD-MET properties (Tables S6 and S7). Of note, selected bioflavonoids were found to be non-inhibitor of several cytochromes (CY) such as CYP2D6, CYP1A2, CYP2C19, CYP2C9, CYP2D6, and CYP3A4, which plays a crucial role in the drug metabolism as well as various xenobiotics; an inhibition of these cytochromes may lead to the reduced drug efficacy, drug activation, and drug metabolism. Also, the selected four bioflavonoids exhibit low gastrointestinal absorption along with a lack of Blood-Brain Barrier (BBB) permeability. However, Rutin and Nicotiflorin showed three violations while Isoquercitrin and Hyperoside showed two violations for the Lipinski's rule of (Tables S6 and S7). The selected bioflavonoids also showed violations against several other rules related to drug-likeness, such as Ghose, Veber, Egan, and Muegge. However, the rules for drug-likeness are not mandatory to be fulfilled by natural compounds as cells distinguish the bioactive compounds through the active transport system [90,91]. Additionally, several other properties related to medicinal chemistry and pharmacokinetics were evaluated for potential four compounds (Tables S6 and S7). Importantly, all of the bioflavonoids showed the negative AMES toxicity test and non-carcinogenic profiles via admetSAR server. Moreover, Conclusively, ADMET analysis suggested the selected bioflavonoids against the ZIKV pro and ZIKV RdRp proteins with ideal medicinal properties.

Long Interval Molecular Dynamics Simulation
In the field of drug discovery using computational approaches, MD simulation is an imperative technique by which the dynamic stability and the formation of intermolecular interactions of docked protein-ligand complexes are analyzed with respect to time.
In this study, four common bioflavonoids, viz. Rutin, Nicotiflorin, Isoquercitrin, and Hyperoside, as bispecific inhibitors, showed substantial docking energy with targeted ZIKV proteins, and the resultant complexes; ZIKV pro −Rutin and ZIKV RdRp −Rutin, ZIKV pro −Nicotiflorin and ZIKV RdRp −Nicotiflorin, ZIKV pro − Isoquercitrin and ZIKV RdRp − Isoquercitrin, ZIKV pro −Hyperoside, and ZIKV RdRp −Hyperoside, were considered for 500 ns explicit MD simulation under constant pressure and temperature to analyze their dynamic stability and intermolecular interaction profiles with respect to time. Additionally, ZIKV pro −O7N inhibitor and ZIKV RdRp −Sofosbuvir inhibitor complexes were also studied under similar MD simulation conditions and marked as reference trajectories for comparative analysis with that of docked complexes of viral proteins with bioflavonoids ( Figures 5-7).       At first, stability and steadiness of docked bioflavonoids in the binding pocket of respective target proteins, i.e., ZIKV pro and ZIKV RdRp , were observed at the end of 500 ns simulation by comparison to the respective initial frames, revealed the acceptable change in conformation of protein structure and deviations in docked ligands positions, similar to the respective reference inhibitors (Figures 5 and 6). Additionally, analysis of intermolecular interaction profiles extracted for the last frames of viral protein-bioflavonoids and reference inhibitors were also found to maintain the several conserved molecular contacts by comparison to their respective initial poses (Table 1, Tables S4 and S5). Altogether, these observations suggested that docked bioflavonoids have substantially occupied the binding pocket of respective viral proteins by comparison to the reference inhibitors during the simulation interval. Hence, the generated 500 ns trajectories of the viral proteins docked with selected bioflavonoids were then analyzed for the statistical analysis, including root mean square deviation (RMSD), root mean square fluctuation (RMSF), and total interaction fraction for protein-ligand (PL) contact mapping by comparison to the respective reference trajectories as function of 500 ns interval.

RMSD and RMSF Analysis
Initially, the protein and ligand RMSD values for the docked complexes of potential bioflavonoids with ZIKV pro and ZIKV RdRp were analyzed with respect to the initial pose as a reference frame (Figure 7). In each docked viral protein-bioflavonoids complexes, ZIKV pro and ZIKV RdRp showed considerable deviations (<3 Å) similar to the protein in respective reference trajectories, i.e., ZIKV pro −O7N inhibitor and ZIKV RdRp −Sofosbuvir inhibitor complexes. These observations were also supported by the respective RMSF values (<2.8 Å), except in the C-terminals (>5 Å) of proteins which can be ignored due to far distance from the binding pockets ( Figure S1). Thus, RMSD analysis of viral proteins reveals no substantial effect of docked bioflavonoids on the global minima of the ZIKV proteins during 500 ns MD simulation interval.
Also, calculated RMSD values for all of the docked bioflavonoids in ZIKV pro showed jumps to higher deviations; Rutin (<6 Å), Nicotiflorin (<10 Å), Isoquercitrin (<12 Å), and Hyperoside (<5 Å) by comparison to O7N inhibitor (<6 Å), followed by a state of steadiness during simulation interval. Similarly, higher deviations in the docked bioflavonoids were noticed with the ZIKV RdRp , Rutin (<5 Å), Nicotiflorin (<6 Å), Isoquercitrin (<7 Å), and Hyperoside (<7 Å) by comparison to Sofosbuvir inhibitor (<5 Å), was notice on some intervals followed by state of equilibrium till the end of the simulation. Interestingly, Rutin (<3 Å in ZIKV pro and <4 Å in ZIKV RdRp ) and Isoquercitrin (<3 Å in ZIKV pro and <5 Å in ZIKV RdRp ) were observed with most acceptable RMSD values and state of global minima by comparison to the reference inhibitors (<6 Å for both O7N inhibitor and Sofosbuvir inhibitor for ZIKV pro and ZIKV RdRp , respectively) at the end of the 500 ns MD simulation interval. The observed RMSD results for the bioflavonoids and reference inhibitors of ZIKV proteins were supported by calculated RMSF values (<4 Å) as a function of simulation interval ( Figure S2). Notably, higher deviations in both the docked bioflavonoids and reference inhibitors were suggested due to the interaction of heavy atoms in the ligands with the active residues in the binding pockets of the viral proteins that resulted in the conformational change in the binding poses of docked ligands, as observed in Figures 5 and 6.

Protein-Ligand Interaction Profiling
In the protein-ligand interaction, non-covalent interactions, especially H-bond and other interactions, such as hydrophobic interaction, ionic interactions, π-π stacking, salt bridges, and water bridges formation, have been reported as essential forces to maintain the stability to the complex. Therefore, in addition to the RMSD and RMSF analysis on the 500 ns simulation trajectory of each docked complex, protein-ligand contact profiling based on non-covalent interactions was also measured for all of the ZIKV pro −bioflavonoids and ZIKV RdRp −bioflavonoids complexes, and compared to the respective reference com-plexes, i.e., ZIKV pro −O7N inhibitor and ZIKV RdRp −Sofosbuvir inhibitor, respectively (Figures 8 and 9, and Figures S3-S7). and ZIKV RdRp −bioflavonoids complexes, and compared to the respective reference complexes, i.e., ZIKV pro −O7N inhibitor and ZIKV RdRp −Sofosbuvir inhibitor, respectively (Figures 8 and 9, and Figures S3-S7).    In the case of ZIKV pro −ligand docked complexes, all of the selected bioflavonoids, i.e., Rutin, Nicotiflorin, Isoquercitrin, and Hyperoside, displayed significant intermolecular contact formation against reference O7N inhibitor (Figure 8). In comparison to the initial ZIKV pro −ligand docked complexes where the residues of ZIKV pro (A:Ser 81 (Figure 3), the most of the residues were also noted in the protein-ligand contact maps obtained from 500 ns MD simulation trajectories, which suggested the stability of docked bioflavonoid in the binding pocket of ZIKV pro . Notably, B:His 51 residue was observed for the formation of hydrophobic interaction for more than 50% of total interaction fraction in all ZIKV pro −bioflavonoid complexes except with Rutin (Figure 8). Similarly, B:TYR 161 residue was observed to form hydrophobic interactions for more than 40% in ZIKV pro −Rutin and ZIKV pro −Nicotiflorin complexes while 90% in ZIKV pro −Isoquercitrin and ZIKV pro −Hyperoside complexes of total interaction fraction. Also, B:Gly 153 was observed to form mainly H-bond in the ZIKV pro −Rutin and ZIKV pro −Nicotiflorin complexes only, and B:Tyr 130 appeared as a naïve interacting residue in all of the ZIKV pro −bioflavonoids complexes during the simulation interval. Moreover, other essential residues, such as A:Ser 81  formed single H-bond for more than 50% of the total interaction fraction calculated as function of 500 ns simulation interval. Also, B:Val 155 and B:Tyr 161 residues showed hydrophobic interactions for more than 50% of the 500 ns simulation time ( Figure S3). Additionally, contribution of water bridge formation was also noted in all of the ZIKV pro complexes with the docked bioflavonoids and reference inhibitor (Figure 8 and Figure S3).

Endpoint Free Binding Energy Calculation
Molecular Mechanics Generalized Born Surface Area (MM/GBSA) method was applied to calculate the net binding free energy on the extracted poses from the last 10 ns MD simulation trajectory of respective docked complexes. Besides, energy dissociation components were also calculated to predict their contribution to the net stability of identified potential bioflavonoids complexes with viral proteins, i.e., ZIKV pro and ZIKV RdRp . The analysis of net binding free energy for the screened bioflavonoids docked with the ZIKV pro and ZIKV RdRp showed considerable energy values by comparison to the respective reference complexes. Interestingly, Rutin docked with the ZIKV pro and ZIKV RdRp showed the highest negative free binding energy compared to other identified bioflavonoids (Table S8, Figure 10). Notably, Rutin was also identified for the formation of stable complex from 500 ns MD simulation via strong intermolecular interactions (Figures 7-9).  Figures S3-S7). Collectively, protein−ligand contact mapping suggests the stability of docked complexes, essentially contributed by the formation of H-bonds and hydrophobic interactions during 500 ns MD simulations.

Endpoint Free Binding Energy Calculation
Molecular Mechanics Generalized Born Surface Area (MM/GBSA) method was applied to calculate the net binding free energy on the extracted poses from the last 10 ns MD simulation trajectory of respective docked complexes. Besides, energy dissociation components were also calculated to predict their contribution to the net stability of identified potential bioflavonoids complexes with viral proteins, i.e., ZIKV pro and ZIKV RdRp . The analysis of net binding free energy for the screened bioflavonoids docked with the ZIKV pro and ZIKV RdRp showed considerable energy values by comparison to the respective reference complexes. Interestingly, Rutin docked with the ZIKV pro and ZIKV RdRp showed the highest negative free binding energy compared to other identified bioflavonoids (Table S8, Figure 10). Notably, Rutin was also identified for the formation of stable complex from 500 ns MD simulation via strong intermolecular interactions (Figures 7-9).  By compared to the net binding free energy, i.e., −74.37 ± 6.44 kcal/mol of reference ZIKV pro −O7N complex, all of the ZIKV pro −bioflavonoids docked complexes showed less binding free energy but in the considerable range (Supplementary Table S8 and Figure 10). In contrast, the binding free energy of ZIKV RdRp −Rutin and ZIKV RdRp −Nicotiflorin docked complexes showed higher values, whereas ZIKV RdRp -Isoquercitrin, and ZIKV RdRp -Hyperoside showed less but close to the binding free energy of ZIKV RdRp −Sofosbuvir inhibitor complex (−59.83 ± 3.85 kcal/mol). In addition, the dissociation energy components were also calculated for all of the docked complexes, where ∆G Bind Lipo , and ∆G Bind vdW contributed to the complex stability, whereas ∆G Bind Solv GB was responsible for destabilizing the respective complexes. (Figure 10 and Figure S8 and Supplementary Table S8). Conclusively, these results suggested that the affinity and stability of ZIKV pro were higher for Rutin, followed by Isoquercitrin, Hyperoside, and Nicotiflorin. In contrast, the stability and affinity of ZIKV RdRp were higher with Rutin followed by Nicotiflorin, Hyperoside, and Isoquercitrin. Hence, net binding free energy calculated using MMGBSA analysis supports the screened four identified common bioflavonoids as putative inhibitors of viral proteins, viz. ZIKV pro and ZIKV RdRp .

Conclusions
The essential role of ZIKV pro and ZIKV RdRp in polyprotein processing and genome replication of ZIKV as well as lack of respective homologs in humans, marks these viral proteins as molecular targets for the development of anti-ZIKV therapeutics. Azadirachta indica plant has placed itself in a category of natural resources with the best medicinal values. Therefore, this study evaluated the reported bioflavonoids from Azadirachta indica for their potential and therapeutic activity against the ZIKV pro and ZIKV RdRp domain using molecular docking simulations, drug-likeness, molecular dynamics simulation, and endpoint binding free energy calculations. Recently, the application of flavonoids as promising source of anti-ZIKV compounds were discussed to exert antiviral activity [92]. Notably, all of the selected bioflavonoids, i.e., Rutin, Nicotiflorin, Isoquercitrin, and Hyperoside, among 44 bioactive compounds against the ZIKV pro and ZIKV RdRp as bispecific inhibitors exhibit considerable binding affinity and dynamic stability. The screened compounds occupied the binding pockets via hydrogen and hydrophobic interactions along with π-π interactions with the essential residues of ZIKV pro and ZIKV RdRp against respective reference inhibitors. The analysis from MD simulation concluded that Rutin and Isoquercitrin with minimum deviation was more stable, followed by Hyperoside and Nicotiflorin with both the viral proteins, i.e., ZIKV pro, and ZIKV RdRp . At last, end-point binding free energy calculation supports the Rutin as potent bispecific inhibitor of ZIKV pro, and ZIKV RdRp . Overall, the present suggested the predicted bioflavonoids from Azadirachta indica as hit candidates and further accurate experimental validation is required to assess their potential as bispecific inhibitors of ZIKV pro , and ZIKV RdRp for the treatment of ZIKV infection.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules27082562/s1, Table S1: List of bioflavonoids reported in Azadirachta indica, used for structure based virtual screening against ZIKVpro and ZIKVRdRp, Table S2: List of virtually screened bioflavonoids from Azadirachta indica against ZIKVpro, Table S3: List of virtually screened bioflavonoids from Azadirachta indica against ZIKVRdRp, Table S4: Intermolecular interactions noted for the screened compounds with the viral proteins, i.e., ZIKVpro and ZIKVR-dRp, within 4 Å around the docked ligand in the respective binding pockets, Table S5: List of various interactions and interacting residues in the active pocket of ZIKVpro and ZIKVRdRp with the selected bioflavonoids were logged from the last pose of respective 500 ns MD trajectories, Table S6: ADME profiling for the selected bioflavonoids from Azadirachta indica as inhibitor against ZIKVpro and ZIKVR-dRp obtained from the swissADME online server (http://www.swissadme.ch/), Table S7: ADMET profiling for the selected bioflavonoids from Azadirachta indica as inhibitor against ZIKVpro and ZIKVRdRp obtained from the admetSAR online server (http://lmmd.ecust.edu.cn/admetsar2/), Table S8: Calculated energy components and net binding free energies (kcal/mol) values for ZIKVpro and ZIKVRdRp complex with selected bioflavonoids against reference compound snap-shots collected from the respective 500 ns MD simulation trajectories, Figure S1: RMSF values plotted for alpha carbon atoms of ZIKVpro and ZIKVRdRp docked with selected bioflavonoids i.e., (a,b) Rutin, (c,d) Nicotiflorin, (e,f) Isoquercitrin, (g,h) Hyperoside, and as well as the reference inhibitors (i) O7N (ZIKVpro reference inhibitor) and (j) Sofosbuvir (ZIKVRdRp reference inhibitor), were extracted from the respective 500 ns MD simulation interval, Figure S2: RMSF values plotted for the bioflavonoids in the docked complexes, i.e., (a) ZIKVpro-Rutin, (b) ZIKVpro-Nicotiflorin, (c) ZIKVpro-Isoquercitrin, (d) ZIKVpro-Hyperoside, (e) ZIKVpro-O7N (Control), (f) ZIKVRdRp-Rutin, (g) ZIKVRdRp-Nicotiflorin, (h) ZIKVRdRp-Isoquercitrin, (i) ZIKVRdRp-Hyperoside, (j) ZIKVRdRp-Sofosbuvir (Control), fit with protein extracted from the respective 500 ns MD simulation interval, Figure S3: Protein-ligand interactions mapping for reference docked complexes, i.e., (a,b) ZIKVpro-O7N, (c,d) ZIKVRdRp-Sofosbuvir inhibitor, extracted from 500 ns MD simulations. In 2D interaction diagram, the residues tyrosine, Valine, and Phenylalanine (green), Aspartic acid (red), histidine and asparagine (blue), and glycine (grey) exhibit the hydrophobic, negative, polar, and non-polar interactions, respectively along with hydrogen bonding (pink arrow) and pi-pi stacking (green line) with the receptor are extracted at 30% of the total MD simulation interaction interval, Figure S4: The panel shows which residues of ZIKVpro interaction with the selected ligands i.e., a. Rutin, b. Nicotiflorin, c. Isoquercitrin, and d. Hyperoside, in each trajectory frame, over the course of the 500 ns md trajectory. Some residues make more than one specific contact with the ligand, which is represented by a darker shade of orange, according to the scale to the right of the plot, Figure S5: The panel shows which residues of ZIKVRdRp interact with the selected ligands, i.e., a. Rutin, b. Nicotiflorin, c. Isoquercitrin, and d. Hyperoside, in each trajectory frame, over the course of the 500 ns md trajectory. Some residues make more than one specific contact with the ligand, which is represented by a darker shade of orange, according to the scale to the right of the plot, Figure S6: 2D interaction diagram of protein-ligand interactions maps for ZIKVpro with selected bioflavonoids, i.e., (a) Rutin, (b) Nicotiflorin, (c) Isoquercitrin, and (d) Hyperoside extracted from the total 500 ns MD simulations. Herein, residues tyrosine, Valine, and Phenylalanine (green), Aspartic acid (red), histidine and asparagine (blue), and glycine (grey) exhibit the hydrophobic, negative, polar, and non-polar interactions, respectively, along with hydrogen bonding (pink arrow) and π-π stacking (green line) with the receptor are extracted at 30% of the total MD simulation interaction interval, Figure S7: 2D interaction diagram of protein-ligand interactions mapping for ZIKVRdRp with selected natural compounds, i.e., (a) Rutin, (b) Nicotiflorin, (c) Isoquercitrin, and (d) Hyperoside, extracted from 500 ns MD simulations. Residues tyrosine, Valine, and Phenylalanine (green), Aspartic acid (red), histidine and asparagine (blue), and glycine (grey) exhibit the hydrophobic, negative, polar, and non-polar interactions, respectively, along with hydrogen bonding (pink arrow) and pi-pi stacking (green line) with the receptor are extracted at 30% of the total MD simulation interaction interval, Figure S8: Calculated free energy components and net MM/GBSA binding free energy (kcal/mol) with standard deviation values for extracted snapshots of reference docked complexes, i.e., (a). ZIKVpro-O7N, (b). ZIKVRdRp-Sofosbuvir, from respective 500 ns MD simulation trajectories. Funding: This research is supported by the charitable donation from the Late Sheikh Ibraheem Ahmed Azhar to cover the article processing charge (APC) as a contribution to the scientific research community.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request. fellowship, and to the Director, AIRF, Jawaharlal Nehru University, New Delhi, India for giving access to Schrödinger suite. Also, the authors are thankful for the received technical support from the Institute of Biotechnology of the Czech Academy of Sciences (Institutional Research Concept, RVO: 86652036).