Identifying the Most Potent Dual-Targeting Compound(s) against 3CLprotease and NSP15exonuclease of SARS-CoV-2 from Nigella sativa : Virtual Screening via Physicochemical Properties, Docking and Dynamic Simulation Analysis

: Background: The outbreak of the coronavirus (SARS-CoV-2) has drastically affected the human population and caused enormous economic deprivation. It belongs to the β -coronavirus family and causes various problems such as acute respiratory distress syndrome and has resulted in a global pandemic. Though various medications have been under trial for combating COVID-19, speciﬁc medicine for treating COVID-19 is unavailable. Thus, the current situation urgently requires effective treatment modalities. Nigella sativa , a natural herb with reported antiviral activity and various pharmacological properties, has been selected in the present study to identify a therapeutic possibility for treating COVID-19. Methods: The present work aimed to virtually screen the bioactive compounds of N. sativa based on the physicochemical properties and docking approach against two SARS-CoV-2 enzymes responsible for crucial functions: 3CLpro (Main protease) and NSP15 (Nonstructural protein 15 or exonuclease). However, simulation trajectory analyses for 100 ns were accomplished by using the YASARA STRUCTURE tool based on the AMBER14 force ﬁeld with 400 snapshots every 250 ps. RMSD and RMSF plots were successfully obtained for each target. Results: The results of molecular docking have shown higher binding energy of dithymoquinone (DTQ), a compound of N. sativa against 3CLpro and Nsp15, i.e., − 8.56 kcal/mol and − 8.31 kcal/mol, respectively. Further, the dynamic simulation has shown good stability of DTQ against both the targeted enzymes. In addition, physicochemical evaluation and toxicity assessment also revealed that DTQ obeyed the Lipinski rule and did not have any toxic side effects. Importantly, DTQ was much better in every aspect among the 13 N. sativa compounds and 2 control compounds tested. Conclusions: The results predicted that DTQ is a potent therapeutic molecule that could dual-target both 3CLpro and NSP15 for anti-COVID therapy. atoms. Overall, the results show that DTQ has a higher afﬁnity towards the target proteins compared to all the tested compounds, conﬁrming it as a potential multi-targeting candidate for COVID-19 treatment.


Introduction
Since the beginning of the 21st century, different coronavirus strains have been identified which cause severe pneumonia among humans [1]. COVID-19 is the name specified for an infectious disease instigated by SARS-CoV-2 (2019 novel coronavirus). Moreover, it originated in Wuhan City, China then it metastasized into a global outbreak. This is the first time, post-World War II, that such economic and social distress has been experienced. Altogether, these circumstances compelled the World Health Organization (WHO) to pronounce COVID-19 as a global pandemic [2]. Until now, nearly 198 million COVID-19 cases, alongside 4,224,466 COVID-19 deaths, have been accounted for, and this contagion proliferation is yet undefeatable with no respite in its annihilation trajectory. A significant clinical characteristic of COVID-19 and the primary cause of death in COVID-19 patients is severe acute respiratory distress syndrome (ARDS) [3,4]. As per WHO, it is still a serious issue affecting public health [5].
MERS-CoV and SARS-CoV belong to the β-coronavirus family following 229E, NL63, OC43, HKU1 viruses; SARS-CoV-2 is the seventh known coronavirus. SARS-CoV-2 has an RNA (+) encapsulated inside the lipid membrane. It has a viral genome that is identical to the SARS-CoV, MERS, bat-CoV RatG13 genome, and the bat is considered the natural host of the virus [6][7][8]. The whole genome of the SARS-CoV-2 can be categorized to ORF's (open reading frames) at the 3 end of the genome, coding the structural and accessory proteins. The structural proteins include the Envelope (E), spike glycoprotein (S), membrane (M), and nucleocapsid (N) proteins [9,10]. There are 16 nonstructural proteins (NSP1-NSP16) in the number encoding for enzymes of vital functions such as replication and transcription, including main protease (3CLpro), RNA-dependent RNA polymerase (RdRP), helicase, Exo, and endonucleases (NSP15) which account for the highly conserved regions for targeting the coronavirus [11,12].
Wide varieties of bioactive molecules are available in natural products with potent antiviral efficacy, among which Nigella sativa is envisaged to have antiviral efficacy against coronavirus infection. It was reported to have antiviral efficacy [13,14]. It was also shown to have anti-inflammatory, immunomodulating, and antimicrobial properties [15,16]. N. sativa has potent antiviral activity against HIV-1 replication, herpes simplex and Hepatitis B virus [17], Hepatitis C Virus [18], and influenza virus H5N1 [19] in vitro. The reported antiviral efficacy and medical importance of N. sativa make it a potential therapeutic candidate against SARS-CoV-2.
The volatile oils of N. sativa are linked with its biological activities, which consist of trans-anethole, α-pinene, β-pinene, carvacrol, carvone, p-cymene, dithymoquinone, limonene, longifolene, α-thujene, thymohydroquinone, thymol, and thymoquinone as major active compounds [20][21][22][23]. All these compounds have been reported to have antiviral potential. Trans-anethole has been reported to repress the infectivity of the herpes virus by 90% with a high selectivity index of 160 [24]. The αand β-pinene compounds showed potency against the infectious bronchitis virus via interacting with its nucleocapsid protein [25]. Carvacrol showed marked antiviral activity against the herpes simplex virus 2 by the intracellular RIP3-arbitrated programmed cell necrosis pathway [26], while carvone has shown an immunostimulatory response in adenovirus infected mouse models [27]. A significant (>than 80%) reduction in plaque formation was observed when p-cymene was incubated with the herpes simplex virus 1 [28]. Similarly, longifolene, α-thujene, and thymohydroquinone also showed antiviral potential against polyomavirus, herpesvirus, and zikavirus, respectively [29][30][31]. Thymol and related monoterpenes were also reported to have potential against the herpes simplex virus 1 [32]. However, thymoquinone is considered the most abundant and major compound present in the volatile oil of N. sativa, among all other compounds. Indeed, thymoquinone is known to activate several immunological responses and has shown potent antiviral activity against different viruses such as the Epstein-Barr virus, cytomegalovirus, human immunodeficiency virus, and murine hepatitis C virus [33]. In addition, thymoquinone is present in its dimeric form as well, which is commonly known as dithymoquinone or nigellone. Dithymoquinone has shown progressive effects on humans during renal-hepatic toxicity, diabetes, and cancer [34][35][36].
All these findings prompted us to explore a potential lead from these 13 active compounds of N. sativa through a molecular docking and simulation approach against 3CLpro and NSP15 as a new hope for predicting effective therapeutics for COVID-19 treatment. However, in vitro and in vivo studies are warranted to validate the findings of the present study.

Calculation of Physicochemical Properties and Prediction of Toxicity Potential
The individual phytochemical constituent of N. sativa was assessed for its toxicity and the physicochemical properties using the Osiris DataWarrior property explorer tool ( www.openmolecules.org, accessed on 2 September 2021) [39]. The tool was used to assess the molecular weight, the number of hydrogen bond donors and acceptors, rotatable bonds, and cLogP value. The molecular properties described by these five parameters were used to evaluate the violation of Lipinski's rule [40], which is used to predict the compound pharmacokinetics in the human body and compound suitability for oral formulations. However, the rule is used for proper drug candidate selection and not to predict the pharmacological activity of the compound. Here, cLogP (the logarithm of partition coefficient between n-octanol and water) was calculated for the prediction of hydrophilicity. The lower the cLogP value, the more the compound will be hydrophilic and with better absorption. In the Osiris DataWarrior tool's analysis, the predicted toxicity values were dependent on the comparison of the precalculated investigated molecules with the structures of the tested molecules. Various aspects and effects of the toxicity, including tumorigenicity, mutagenicity, and irritability of the tested compounds, were also studied using the Osiris DataWarrior tool.

Molecular Docking
The molecular docking study was carried out as per the protocol followed by Rizvi et al. [41]. After the addition of gasteiger partial charges to the proteins, the MMFF94 force field was applied for the energy minimization of each compound. The rotatable bonds were defined following the inclusion of non-polar hydrogen atoms. Then the solvation parameters were added with Kollman united atom type charges and hydrogen atoms. The grids were created using the auto grid in the dimension 60 × 60 × 60 Å with the spatial distance between the points of approximately 0.375 Å. For specific targeting, the X, Y, and Z coordinates were fed as −16.539, 15.246, and 67.334 for 3CLpro and −72.419, 24.902, and −37.184 for Nsp 15, respectively. To calculate the electrostatic terms and Van der Waals forces, they were maintained with default parameters and dielectric functions in AutoDock. The molecular docking simulation was conducted using the 'Solis and Wets local search method' and 'Lamarckian genetic algorithm', where other parameters including torsion, orientation, and initial positions were arbitrarily set. Furthermore, 100 different runs were set, which ended after 2,500,000 evaluations for each experiment. The size of the population for the study was 150.

LIGPLOT+ Analysis
The best-docked conformation was chosen and analyzed using LIGPLOT+ Version v.2.1 [42]. The hydrophobic and hydrogen bond interactions were examined and the LIGPLOT algorithm was applied to convert 3D into 2D final figures.

Molecular Dynamics Simulation Study
The docking complexes of 3CLpro and Nsp15 bounded with dithymoquinone (DTQ) (the ligand identified and selected by screening) were subjected to molecular dynamics simulation analysis using YASARA STRUCTURE v.20.12.24.W.64 [43]. The AMBER14 force field was used to set the simulation cell boundary to periodic. The simulation cell was permitted to contain 20 Å surrounding the protein and filled with solvent (water) at a density of 0.997 g/mL. To enhance the stability of the solute, the H bonding network was optimized [44], followed by fine-tuning the target protein protonated state via pKa prediction at 7.4 pH. However, to neutralize the environment, NaCl ions were also added to the system. At first, the energy of the system was minimized in the YASARA structure tool to rectify the covalent geometry and bumps. Here, simulated annealing and steepest descent minimization techniques were applied to eliminate clashes before running the simulation trajectory for 100ns by the 'AMBER14 force field' for the solute [45], 'TIP3P' for water, and 'AM1BCC [46] and GAFF2 [47]' for DTQ. 'Van der Waals forces' were applied with an 8 Å cut-off, while 'electrostatic forces and the Particle Mesh Ewald algorithm' were applied without any cut-off limit. Using NPT ensembles [48], the systems were equilibrated by applying position restrains. Equilibration for both the complexes was integrated with a multiple time step of 2.5 fs and 5.0 fs for bonded and non-bonded interactions, respectively, at 1 bar pressure and a temperature of 298 K. NPT ensembles of isobaric and isothermal conditions were maintained for the entire simulation process. The temperature control rescaling approach was used to initiate a Berendsen thermostat [49], which is based on the time-averaged temperature. Berendsen barostat [49], based on time average pressure, was used to control pressure. However, the tuned version of LINCS [50] via the multiple time-step algorithm [48] was used to constrain bonds and angles involving hydrogen. It is noteworthy to mention that the easy interphase provided by the macros tool in the YASARA structure, i.e., 'md_runfast.mcr' was used to run the total simulation process in the present study. In addition, 'md_analyze.mcr' was used to analyze the trajectory path of the simulation process. Two separate simulations were performed for 'DTQ-3CLpro' and 'DTQ-Nsp15', respectively. Figures for each simulation were generated by using the YASARA structure tool where 400 snapshots were taken each 250 ps. A ray-traced diagram of each simulation system with a cubic shape box was also generated via the YASARA structure tool (Figures S1 and S2).

Virtual Screening
Initially, fifteen compounds including thirteen active components of N. sativa and two controls (lopinavir and benzopurpurin B) were subjected to virtual screening via molecular docking, physicochemical and toxicity properties analysis, followed by the simulation of the selected lead compound (Figure 1). DTQ has shown an excellent binding affinity toward both the enzymes compared to all the other compounds; whereas, the standard lopinavir exhibited a binding affinity of −7.95 kcal/mol against 3CLpro main protease and benzopurpurin B showed a binding affinity of approximately −5.87 kcal/mol against Nsp15. Table 1 shows the binding energy scores for ligands' molecular docking against the viral enzyme proteins 3CLpro and Nsp15.
As per the results, DTQ has been selected as the molecule of interest for effectively combating the SARS-CoV-2. The LigPlot analysis of 'DTQ-3CLPro' interaction in Figure 2a shows DTQ formed hydrogen at the Glu 166 from the O1 molecule and has a bond length of approximately 3.15 Å. Figure 2b shows the 2D binding scheme of the DTQ and Nsp15 interaction, in which two hydrogen bonds are formed between DTQ and the active site of Nsp15 at Ile 169 and Asn200 amino acid residue with a bond length of 3.23 Å and 2.61 Å, respectively. The hydrogen bond was formed between DTQ-Ile169 at O4 and N atoms, and DTQ-Asn200 between O3 and ND2 atoms. Overall, the results show that DTQ has a higher affinity towards the target proteins compared to all the tested compounds, confirming it as a potential multi-targeting candidate for COVID-19 treatment.

Virtual Screening via Physicochemical Properties Analysis
The physicochemical properties and toxicity potential of phytoconstituents of N. sativa were assessed using the Osiris DataWarrior property explorer tool for drug-likeness. The 13 phytoconstituents of N. sativa, including trans-anethole, α-pinene, β-pinene, car-vacrol, carvone, p-cymene, dithymoquinone, limonene, longifoline, α-thujene, thymohydroquinone, thymol, and thymoquinone have not violated the Lipinski rule as shown in Table 2. Among the control drugs, both lopinavir and benzopurpurin violated the Lipinski rule with a molecular weight of approximately 628.81 g/mol and 680.76 g/mol higher than 500 g/mol, respectively. Lopinavir and benzopurpurin also have a higher number of rotatable bonds and hydrogen bond receptors. Our molecule of interest, DTQ, has a molecular weight of 328.40 g/mol, which is suitable for crossing the biosystem's membrane barriers. It also has a cLogP value of approximately 2.73, indicating the compound's moderately lipophilic nature with proper absorption in the gastro-intestinal tract with no violation of the Lipinski rule. Thus, it makes a promising candidate for combating COVID-19. * Control drugs/compounds; ** Logarithm of compound partition coefficient between n-octanol and water.

Virtual Screening via Toxicity Assessment
The toxicity assessment of the phytoconstituents was carried out using the Osiris DataWarrior property explorer tool. Of thirteen compounds, three compounds, including βpinene, dithymoquinone, and longifoline, had toxic effects such as mutagenic, tumorigenic, and reproductive irritant effects (Table 3). Toxicity analysis showed that DTQ remained inactive for the predictive toxicity parameters used and was expected to be non-toxic.

Molecular Dynamic Simulation Study
Based on virtual screening, dithymoquinone (DTQ) was selected for dynamic simulation analysis. Molecular dynamic (MD) simulation studies were performed on the YASARA structure tool v.20.12.24.W.64 (using AMBER14 force field) to understand the protein-ligand complex's dynamic activity in a solvent environment over a period of time. Simulated conditions were optimized to map the stabilization of the 'DTQ-3CLpro' and 'DTQ-NP15' complexes and all the parameters were maintained throughout the simulation time. Simulations were performed three times for a time period of 100ns with the naïve protein alone and with the ligand-protein complex. Understanding the binding affinity and the protein-ligand complex's stability in a time-bound environment was the ultimate aim for performing the dynamic studies. The AMBER14 force field was employed to generate a plot for the total potential energy of the system (Figure 3). Once the simulation is initiated from a ground zero or frozen energy-minimized confirmation, a steep energy rise is typically observed through the initial picoseconds of the simulation run. This rise in energy is due to partial storage of kinetic energy as potential energy and, importantly, the potential energy usually does not decrease during a large time-scale due to counter ions. They are positioned primarily with lower potential energy near the charged solute groups. From here, they will separate to increase potential energy and entropy. It can be observed from the plot (Figure 3a) that the potential energy of the system for the DTQ-3CLpro complex fluctuated from −1,223,500 kJ/mol to −1,218,500 kJ/mol. While, for the DTQ-NSP15 complex (Figure 3b), the total potential energy of the system fluctuated from −1,219,000 kJ/mol to −1,213,000 kJ/mol. Importantly, fluctuations in both the systems were within an acceptable range demonstrating the stability and validity of the simulation.
The solute root mean square deviation (RMSD) for the starting structure was plotted for both the 'DTQ-3CLpro' and 'DTQ-NSP15' complexes ( Figure 4). It shows three RMSDs, i.e., RMSD for Calpha (Ca), RMSD for backbone (Bb), and RMSD for all heavy atoms (All) for each complex. The RMSD results (Figure 4a) for the DTQ-3CLpro complex showed that RMSD Ca and RMSDBb plots have slight fluctuations until 20 ns; although, from 20 ns to 100 ns both the plots were overlaying with minimal fluctuations. However, backbone RMSD fluctuations for the 'DTQ-3CLpro' complex were within the range of between 1 and 2.8 Å.  On the other hand, the RMSD results (Figure 4b) of the 'DTQ-NSP15' complex showed the overlapping of the RMSDCa and RMSDBb plots from the beginning of the simulation until 100 ns. After 65 ns, little fluctuation was observed, while the overall range of backbone RMSD fluctuations was acceptable, i.e., 1.5-2.8 Å. For the 'DTQ-3CLpro' and 'DTQ-NSP15' complexes, RMSD for heavy atoms showed fluctuations between 20 ns and 12 ns within the limit of 9 Å and 8 Å, respectively. Backbone RMSD results depicted that the 'DTQ-3CLpro' complex would reach equilibrium after 35 ns time, whereas the 'DTQ-NSP15' complex would reach equilibrium after 5 ns time (with little fluctuations after 65 ns), reflecting structural protein flexibility being reserved in their respective complex form.
Residual-wise fluctuations could be seen from the root mean square fluctuation (RMSF) plot ( Figure 5), where RMSF/solute residue was evaluated from the (avg.) RMSF of each atom was included in the residue. The RMSF for the DTQ-3CLpro complex (Figure 5a) and DTQ-NSP15 complex (Figure 5b) were found to be 5.18 Å and 3.98 Å, respectively. The radius of gyration for both the simulated protein complexes was plotted in Figure 6. The calculation for the radius of gyration is dependent on the center of mass of the protein that indicates the compactness of the protein structure. It will remain steady if the protein shows stability; however, due to instability, it will change with time. Importantly, in our study, both 3CLpro ( Figure 6a) and Nsp15 (Figure 6b) showed little fluctuation that corresponded to their stability with time. Fluctuations for 3CLpro and Nsp15 were in an acceptable range between 1.6 Å and 1.4 Å, respectively.
Stability and successful protein folding were also analyzed by the number of hydrogen bonds. The hydrogen bond analysis showed a growing number of bonds with respect to time during the entire 100 ns simulation for both simulations (Figure 7). It was apparent that, although the number of hydrogen bonds fluctuated during simulation, the stable nature of the protein did not change.
In addition, the dynamic cross-correlation matrix for both complexes was also obtained from the simulation analysis (data not shown for concision). However, some of the snapshots from the simulation process for the DTQ-3CLpro complex and DTQ-NSP15 complex are shown in Figures 8 and 9, respectively. Conclusively, all the simulation data suggested that both complexes were stable after initial fluctuations; however, DTQ-NSP 15 showed more stability than the DTQ-3CLPro complex.

Discussion
The outbreak of COVID-19 in Wuhan, China's Hubei province, opens up new perspectives for discovering synthetic and natural small compounds to find leads that could be used to develop feasible antiretroviral medications. The sequenced 3D-crystallographic structures of 3CLpro and NSP15 enzyme proteins of SARS-CoV-2 were available from the PDB database [37,38].
The 3CLpro has been responsible for regulating protease activity in the virus, making it a therapeutic target for COVID-19. The protease activity was attributed to the Cys-His catalytic active site, which subsequently leads to the inhibition of the protease activity responsible for replicating the virus [51]. The secondary structure of the 3CLpro protein consists of three domains where the antiparallel β-barrel can be observed in Domains I and II; whereas, in domain III, the antiparallel globular cluster was connected via long loop to Domain II [52,53]. Recently NSP15 has been identified as a therapeutic target for COVID-19 treatment as they are essential for their life cycle and virulence. Each unit of NSP15 is composed of approximately 345 amino acids, which are folded into three domains. They are crucial for virus replicase and transcriptase formation. NSP15 tends to suppress the innate immunity (type I IFN (IFN-α/β)) via infecting macrophages and causing the evasion of the virus from the immune system [54,55] that makes it a susceptible target for COVID-19 therapeutics. The work aimed to identify the potential therapeutic lead from the N. sativa phytoconstituents due to its previously reported pharmacological properties for antiviral treatment. The phytoconstituents were screened by physicochemical properties and molecular docking, followed by dynamics simulation against the therapeutic targets, i.e., 3CLpro and NSP15.
A carbonyl polymer of thymoquinone, dithymoquinone (DTQ), was a more promising compound against the selected targets. The results clearly showed that the DTQ had exerted higher binding affinity towards the two target proteins with a binding energy of −8.56 kcal/mol and −8.31 kcal/mol against 3CLpro and Nsp15, respectively. The binding energy of all phytoconstituents against target proteins is shown in Table 1. The physicochemical properties were assessed using the Osiris DataWarrior property explorer tool. Results from Tables 2 and 3 show that DTQ does not have any toxic effects and obeys the Lipinski rule, which makes it a safe therapeutic candidate. Its partial water solubility character makes it suitable for oral administration. Importantly, DTQ was better than the positive control Lopinavir (3CLpro) and benzopurpurin B (Nsp15) in every aspect. Molecular dynamic simulation for 100 ns with target proteins in their naïve and ligand complex form revealed good stability of DTQ-target proteins that are evident from the RMSD and hydrogen bond formation. Additionally, DTQ has shown that it can target both proteins (3CLpro and Nsp15) with almost equal efficacy, which is predicted to increase its antiviral efficacy against SARS-CoV-2.

Conclusions
The results predict that dithymoquinone (DTQ) has good affinity and stability against the target proteins 3CLpro and NSP15. It can exert its antiviral effect by inhibiting the virus's vital functions such as replication, protease activity, and life cycle. DTQ might be used as a potential therapeutic molecule against COVID-19. However, further in vitro and in vivo evaluation is needed to assess its efficacy for optimal transformation to clinical administration.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/pr9101814/s1, Figure S1: A ray-traced picture of the simulated system for 'DTQ-3CLpro' complex. Here, the simulation cell bound-ary is set to periodic. Figure S2: A ray-traced picture of the simulated system for 'DTQ-Nsp15' complex. Here, the simulation cell boundary is set to periodic.