The Inhibitory Potential of 2′-dihalo Ribonucleotides against HCV: Molecular Docking, Molecular Simulations, MM-BPSA, and DFT Studies

Sofosbuvir is the first approved direct-acting antiviral (DAA) agent that inhibits the HCV NS5B polymerase, resulting in chain termination. The molecular models of the 2′-dihalo ribonucleotides used were based on experimental biological studies of HCV polymerase inhibitors. They were modeled within HCV GT1a and GT1b to understand the structure–activity relationship (SAR) and the binding interaction of the halogen atoms at the active site of NS5B polymerase using different computational approaches. The outputs of the molecular docking studies indicated the correct binding mode of the tested compounds against the active sites in target receptors, exhibiting good binding free energies. Interestingly, the change in the substitution at the ribose sugar was found to produce a mild effect on the binding mode. In detail, increasing the hydrophobicity of the substituted moieties resulted in a better binding affinity. Furthermore, in silico ADMET investigation implied the general drug likeness of the examined derivatives. Specifically, good oral absorptions, no BBB penetration, and no CYP4502D6 inhibitions were expected. Likely, the in silico toxicity studies against several animal models showed no carcinogenicity and high predicted TD50 values. The DFT studies exhibited a bioisosteric effect between the substituents at the 2′-position and the possible steric clash between 2′-substituted nucleoside analogs and the active site in the target enzyme. Finally, compound 6 was subjected to several molecular dynamics (MD) simulations and MM-PBSA studies to examine the protein-ligand dynamic and energetic stability.


Introduction
Hepatitis C virus (HCV) is one of the major causes of chronic liver diseases and hepatocellular carcinoma (HCC) is the third leading cause of cancer death worldwide. The HCV nonstructural protein 5B (NS5B) is an RNA-dependent RNA polymerase (RdRp) that plays a vital role in RNA viral replication. Nucleoside/tide analogs have been used to treat viral infections and cancer for decades. Since the first FDA-approved nucleoside analog "edoxudine" in 1969, there are currently 25 FDA-approved nucleoside/tide analogs and another 15 nucleoside analogs used as antiviral and anticancer agents, respectively. Sofosbuvir (phosphoramidate prodrug of 2 -deoxy-2 -α-fluoro-β-C-methyluridine), a revolutionary HCV NS5B polymerase inhibitor and oral direct-acting antiviral (DAA), one of the best-selling drugs in the world, was approved by FDA in December 2013. Since its FDA approval, the cure rate has risen sharply to 90% after 12 weeks of treatment. Subsequently, the FDA approved many HCV drugs and currently play a pivotal role in HCV cure.
Sofosbuvir is the Sp isomer and has better activity than the Rp isomer [7]. Oral sofosbuvir is absorbed rapidly (Tmax ≈ 1 h), and all inactive metabolites are cleared mainly by the kidney. Sofosbuvir is phosphorylated to its active triphosphate in the hepatocytes, firstly by hydrolysis of the carbonyl ester in the phosphoramidate moiety by the hepatically expressed carbonyl esterase 1 (hCE1) and Cathepsin A (CatA). Further hydrolysis by histidine triad nucleotide-binding protein 1 (H1NT1) to its parent nucleoside (the monophosphate form) is further phosphorylated by phosphorylation kinases to nucleoside diphosphate and then to the triphosphate. This active form can compete with endogenous nucleotides, leading to the termination of HCV-RNA replication ( Figure 2) [7][8][9]. Prodrugs are an inactive or less active form of pharmaceutical compounds, which are metabolized into an active form after administration to improve the delivery to the cell. Pronucleotides are prodrugs in which two different chemical linkers mask the monophosphate moiety, either oxygen to form phosphoester (phosphorous-oxygen bond) or nitrogen to form phosphoramidate (phosphorous-nitrogen bond, The ProTide Technology) prodrugs; phosphorous-nitrogen bonds in phosphoramidate prodrugs are intracellularly cleaved by phosphoramidase [4][5][6]. In 2013, the FDA approved the first direct-acting antiviral drug, Sovaldi ® (sofosbuvir), remarkably treating and curing HCV infection. Sofosbuvir is a phosphoramidate prodrug inhibitor of HCV polymerase.
Sofosbuvir is the Sp isomer and has better activity than the Rp isomer [7]. Oral sofosbuvir is absorbed rapidly (T max ≈ 1 h), and all inactive metabolites are cleared mainly by the kidney. Sofosbuvir is phosphorylated to its active triphosphate in the hepatocytes, firstly by hydrolysis of the carbonyl ester in the phosphoramidate moiety by the hepatically expressed carbonyl esterase 1 (hCE1) and Cathepsin A (CatA). Further hydrolysis by histidine triad nucleotide-binding protein 1 (H1NT1) to its parent nucleoside (the monophosphate form) is further phosphorylated by phosphorylation kinases to nucleoside diphosphate and then to the triphosphate. This active form can compete with endogenous nucleotides, leading to the termination of HCV-RNA replication ( Figure 2) [7][8][9].
Halogen atoms have a remarkable mimicking effect on many atoms and groups. The incorporation of one or more halogen atoms, especially the fluorine atom, into nucleoside analogs improved their pharmacological properties [10,11]. Recently, multiple research teams have discovered and developed novel 2 -deoxy-2 ,2 -dihalonucleosides/tides with high potency against HCV NS5B polymerase [12][13][14][15][16][17][18]. Some of those analogs (e.g., compounds 6, 8, and 12) displayed better potency as compared to the marketed drug Sofosbuvir [12] due to the similarity of the van der Waal's radius of the chlorine atom or bromine with the methyl group, which could provide the steric interference and the polarizability to terminate the RNA chain elongation of the viral replication (Table 1). This work aimed to study and understand the binding modes between the active 2 -deoxy-2 ,2 -dihalonucleotides as a ligand, and NS5B polymerase as a receptor using modern computational techniques to realize the structural-activity relationship and to predict the Table 1. HCV GT1a and GT1b replicon activity and cytotoxicity for the prodrugs 1-13 [12]. copyright has been obtained from Elsevier BMC, License Number 5346420180398, and it will be provided.

Docking Studies
Molecular docking studies were performed to obtain further insights into the binding modes and orientations of the examined compounds in the HCV RNA-dependent RNA polymerase or nonstructural protein 5B (NS5B), specifically the GT-1a and GT-1b subtypes (HCV NS5B GT1a and HCV NS5B GT1b) [19][20][21]. The target proteins were retrieved from Protein Data Bank (PDB ID for GT1a: 4khm and PDB ID for GT1b: 4kai). The co-crystalized ligand (IPV) was used as a reference ligand. This study evaluated how our small molecules (ligands) and the target macromolecules (HCV NS5B GT1a and HCV NS5B GT1b) fit together. Docking studies were carried out using MOE Software. The structures of the tested compounds and the co-crystallized ligand are presented in Figure 3. The docking scores are summarized in Table 2. From the results presented in Table 1, compounds 6   field and ASE as a scoring function. The validation step proved the suitability of the used protocol for the intended docking studies, as demonstrated by the small RMSD (1.12 and 1.23 Å for GT1a and GTIb, respectively) between the docked poses and the co-crystallized ligands. Moreover, the ability of the docking algorithm to retrieve the reported binding mode of the co-crystallized ligands also confirmed the validity of the selected docking algorithm (Figure 4).

Docking Protocol Validation
Validation of the docking protocol was achieved through re-docking of the co-crystallized ligands in the HCV NS5B GT1a and HCV NS5B GT1b active sites using MMFF94X as a force field and ASE as a scoring function. The validation step proved the suitability of the used protocol for the intended docking studies, as demonstrated by the small RMSD (1.12 and 1.23 Å for GT1a and GTIb, respectively) between the docked poses and the cocrystallized ligands. Moreover, the ability of the docking algorithm to retrieve the reported binding mode of the co-crystallized ligands also confirmed the validity of the selected docking algorithm ( Figure 4).

Docking Studies against GT1a
Firstly, the reported data indicated that six amino acids constitute the active site of HCV-NS5B GT1a, and they were involved in the interactions with HCV-NS5B inhibitors. These interacting residues are Arg200, Ser365, Gly192, Tyr316, Phe193, and Tyr448 [22].
The co-crystalized ligand (IPV) displayed a −9.16 kcal/mol binding energy value against GT1a. The detailed binding mode was as follows: the (2-fluorophenyl)boronic acid moiety formed two hydrophobic interactions with Met414 and Cys366. Moreover, the fluoro atom formed a hydrogen bond with the amino acid Tyr448. Furthermore, the cyclohexyl moiety formed three hydrophobic interactions with Cys316, Phe193, and Tyr448. Additionally, the methylsulfonamide group formed two hydrogen bonds with Arg200. The N-methylbenzofuran-3-carboxamide moiety was incorporated into two hydrophobic interactions with Cys336 and Leu384. Moreover, it formed two hydrogen bonds with Arg200 and Ser365. Finally, the terminal 4-fluorophenyl moiety formed hydrophobic stacking with Ile363, Leu204, and Val321 ( Figure 5).

Docking Studies against GT1a
Firstly, the reported data indicated that six amino acids constitute the active site of HCV-NS5B GT1a, and they were involved in the interactions with HCV-NS5B inhibitors. These interacting residues are Arg200, Ser365, Gly192, Tyr316, Phe193, and Tyr448 [22].
The co-crystalized ligand (IPV) displayed a −9.16 kcal/mol binding energy value against GT1a. The detailed binding mode was as follows: the (2-fluorophenyl)boronic acid moiety formed two hydrophobic interactions with Met414 and Cys366. Moreover, the fluoro atom formed a hydrogen bond with the amino acid Tyr448. Furthermore, the cyclohexyl moiety formed three hydrophobic interactions with Cys316, Phe193, and Tyr448. Additionally, the methylsulfonamide group formed two hydrogen bonds with Arg200. The N-methylbenzofuran-3-carboxamide moiety was incorporated into two hydrophobic interactions with Cys336 and Leu384. Moreover, it formed two hydrogen bonds with Arg200 and Ser365. Finally, the terminal 4-fluorophenyl moiety formed hydrophobic stacking with Ile363, Leu204, and Val321 ( Figure 5). Compound 1 (sofosbuvir) exhibited a −7.19 kcal/mol binding affinity value against GT1a. The triphosphate moiety formed five hydrogen bonds with Arg200, Cys336, and Met414. Additionally, it formed a hydrophobic interaction with Phe415. The hydroxyl group of the 2′-deoxy-2′-α-fluoro-2′-β-methylribose moiety formed two hydrogen bonds with the amino acids Ser368 and Ser365. Moreover, it formed two hydrophobic interactions with Ile363 and Val321. The uracil base formed two hydrogen bonds with Leu314 and Cys316. Furthermore, it formed four hydrophobic interactions with Cys366, Cys316, Val321, and Leu314 ( Figure 6). Compound 1 (sofosbuvir) exhibited a −7.19 kcal/mol binding affinity value against GT1a. The triphosphate moiety formed five hydrogen bonds with Arg200, Cys336, and Met414. Additionally, it formed a hydrophobic interaction with Phe415. The hydroxyl group of the 2 -deoxy-2 -α-fluoro-2 -β-methylribose moiety formed two hydrogen bonds with the amino acids Ser368 and Ser365. Moreover, it formed two hydrophobic interactions with Ile363 and Val321. The uracil base formed two hydrogen bonds with Leu314 and Cys316. Furthermore, it formed four hydrophobic interactions with Cys366, Cys316, Val321, and Leu314 ( Figure 6). The binding mode of the active form of compound 6 is illustrated in Figure 7. It showed a binding energy of −6.21 kcal/mol against GT1a. The triphosphate moiety formed three hydrogen bonds with Ser368 and Arg200 in addition to electrostatic interaction with Phe415. The uracil base formed four hydrogen bonds with Arg394, Glu143, and Lys141. The binding mode of the active form of compound 6 is illustrated in Figure 7. It showed a binding energy of −6.21 kcal/mol against GT1a. The triphosphate moiety formed three hydrogen bonds with Ser368 and Arg200 in addition to electrostatic interaction with Phe415. The uracil base formed four hydrogen bonds with Arg394, Glu143, and Lys141. Compound 8 showed a binding energy of −6.79 kcal/mol against GT1a. The triphosphate moiety formed five hydrogen bonds with Asp319, Cys316, and Arg200 and two electrostatic interactions with Asp319. The fluoro atom of the 2′-deoxy-2′-β-chloro-2′-αfluororibose moiety formed a hydrogen bond with Tyr448. Additionally, the uracil base formed another hydrogen bond with Arg386 ( Figure 8). Compound 8 showed a binding energy of −6.79 kcal/mol against GT1a. The triphosphate moiety formed five hydrogen bonds with Asp319, Cys316, and Arg200 and two electrostatic interactions with Asp319. The fluoro atom of the 2 -deoxy-2 -β-chloro-2 -αfluororibose moiety formed a hydrogen bond with Tyr448. Additionally, the uracil base formed another hydrogen bond with Arg386 ( Figure 8).
The co-crystalized ligand (IPV) showed a −7.24 kcal/mol binding energy value against GT1b. Each (2-fluorophenyl) boronic acid moiety and cyclohexyl moiety formed one hydrophobic interaction with Cys366 and Tyr448, respectively. The fluoro atom formed a hydrogen bond with Ser552. Additionally, the N-methyl benzofuran-3-carboxamide moiety was incorporated in four hydrophobic interactions with Cys366 and Leu384 in addition to forming a hydrogen bond with Arg200. Finally, the terminal 4-fluorophenyl moiety formed hydrophobic stacking with Ile363, Leu204, and Val321( Figure 9).
The co-crystalized ligand (IPV) showed a −7.24 kcal/mol binding energy value against GT1b. Each (2-fluorophenyl) boronic acid moiety and cyclohexyl moiety formed one hydrophobic interaction with Cys366 and Tyr448, respectively. The fluoro atom formed a hydrogen bond with Ser552. Additionally, the N-methyl benzofuran-3-carboxamide moiety was incorporated in four hydrophobic interactions with Cys366 and Leu384 in addition to forming a hydrogen bond with Arg200. Finally, the terminal 4-fluorophenyl moiety formed hydrophobic stacking with Ile363, Leu204, and Val321( Figure 9). Compound 1 (sofosbuvir) exhibited a −6.30 kcal/mol binding affinity value against GT1b. The triphosphate moiety formed three hydrogen bonds with Cys366 and Asp319. Additionally, it formed two hydrophobic interactions with Asp319. The hydroxyl group and fluoro atom of the 2′-deoxy-2′-α-fluoro-2′-β-methylribose moiety formed two hydrogen bonds with Asn316 and Arg200. In addition, it formed two hydrophobic interactions with Leu314 and Val321. The uracil base formed two hydrophobic interactions with Val321 and Ile363 ( Figure 10). Compound 6 showed a binding energy of −6.63 kcal/mol against GT1b. The triphosphate moiety formed four hydrogen bonds with Arg200, Try448, Ser368, and Asn316, in addition to electrostatic interaction with Tyr415. The 2′-deoxy-2′,2′-dichlororibose moiety was incorporated in two hydrophobic interactions with Met414 and Tyr415. Finally, the uracil base formed two hydrogen bonds with Arg158 and Arg386 ( Figure 11). Compound 6 showed a binding energy of −6.63 kcal/mol against GT1b. The triphosphate moiety formed four hydrogen bonds with Arg200, Try448, Ser368, and Asn316, in addition to electrostatic interaction with Tyr415. The 2 -deoxy-2 ,2 -dichlororibose moiety was incorporated in two hydrophobic interactions with Met414 and Tyr415. Finally, the uracil base formed two hydrogen bonds with Arg158 and Arg386 ( Figure 11 Compound 8 showed a binding energy of −6.39 kcal/mol against GT1b. The triphosphate moiety formed seven hydrogen bonds with Arg386, Arg158, Asp319, and Cys366, along with an electrostatic attraction with Asp319. The chloro atom of the 2′-deoxy-2′-βchloro-2′-α-fluororibose moiety formed two hydrophobic interactions with Met414 and Tyr448. Furthermore, the uracil base formed two hydrophobic interactions with Cys366 and Tyr415 ( Figure 12). In conclusion, all examined compounds exhibited a good binding mode against the target receptors. It was noticed that the change in the substitution at the ribose sugar produced a mild effect on the binding potential. Specifically, it was found that the fluoro atom has more advantages than the other corresponding groups (chloro, bromo, and methyl). These advantages were confirmed by the increased biological activity of compounds substituted with a fluoro atom. In addition, as appeared in compound 8 (Figure 8), the fluoro atom has the ability to form a hydrogen bond at the active site, producing a significant fitting. Moreover, the chloro atom showed a good fitting with the receptors via its hydrophobic interactions as shown in compound 8 ( Figure 12). From the binding mode of compound 1 (sofosbuvir) (Figures 6 and 10), it is obvious that the methyl group has a significant role in the binding with the receptor via the formation of many hydrophobic interactions. In conclusion, all examined compounds exhibited a good binding mode against the target receptors. It was noticed that the change in the substitution at the ribose sugar produced a mild effect on the binding potential. Specifically, it was found that the fluoro atom has more advantages than the other corresponding groups (chloro, bromo, and methyl). These advantages were confirmed by the increased biological activity of compounds substituted with a fluoro atom. In addition, as appeared in compound 8 (Figure 8), the fluoro atom has the ability to form a hydrogen bond at the active site, producing a significant fitting. Moreover, the chloro atom showed a good fitting with the receptors via its hydrophobic interactions as shown in compound 8 ( Figure 12). From the binding mode of compound 1 (sofosbuvir) (Figures 6 and 10), it is obvious that the methyl group has a significant role in the binding with the receptor via the formation of many hydrophobic interactions.

In Silico ADMET Analysis
In silico ADMET studies were carried out using sofosbuvir as a reference drug. The designed ADMET studies include the following descriptors: (i) blood-brain barrier penetration, which predicts the ability of a molecule to penetrate the blood-brain barrier; (ii) intestinal absorption, which predicts human intestinal absorption of a compound (HIA) after oral administration; (iii) aqueous solubility, which predicts the solubility of a compound in water at 25 • C; (iv) CYP2D6 binding predicts a molecule's potential to inhibit the cytochrome P450 2D6 enzyme in the human body; and (v) plasma protein binding computes the fraction ratio of a molecule that is bound to plasma proteins in the bloodstream [23]. Discovery studio 4.0 software was utilized to investigate the ADMET descriptors for the presented compounds. The predicted ADMET descriptors are documented in Table 3. The results revealed that all compounds exhibited a low or very low BBB penetration ability. Accordingly, such compounds were expected to be safe to CNS. All compounds were expected to have an optimal range of ADMET aqueous solubility levels. Intestinal absorption is the percentage of the absorbed compound from the human gut wall [24]. A well-absorbed compound is absorbed with a level of 90% or more into the bloodstream [25]. According to the in silico ADMET studies, all compounds were predicted to have good intestinal absorption levels.
The cytochrome P450 2D6 (CYP2D6) model computes the potential of a drug to inhibit the CYP2D6 enzyme by utilizing the 2D chemical structure of that drug as input. CYP2D6 is a vital enzyme responsible for the metabolism of a vast range of compounds inside the liver. Accordingly, its inhibition will cause a serious drug-drug interaction. Hence, the expectation of the CYP2D6 inhibition potential is demanded in the field of drug discovery and development [26]. All compounds were anticipated to be non-inhibitors of CYP2D6. Hence, no liver dysfunction side effect is predicted after its administration. The plasma protein-binding model calculates whether a molecule will bind highly (>=90%) to carrier plasma proteins in the blood or not [27]. The output of this experiment could classify whether a molecule is highly bound (>=90% bound) or not to plasma proteins by utilizing the Bayesian cutoff score of −2.209. All compounds were expected to bind plasma proteins with a level of less than 90% (Figure 13). sify whether a molecule is highly bound (>=90% bound) or not to plasma proteins by utilizing the Bayesian cutoff score of −2.209. All compounds were expected to bind plasma proteins with a level of less than 90% ( Figure 13).

Figure 13.
ADMET_AlogP98: lipid-water partition coefficient; ADMET_PSA_2D: polar molecular surface area. Each drug is plotted with a 2D polar surface area (PSA_2D) against its computed atomtype partition coefficient (ALogP98). The area encompassed by the ellipse represents good absorption without any violation of the ADMET properties. The 95% and 99% confidence limit ellipses corresponding to the blood-brain barrier (BBB) and intestinal absorption models are indicated.

Toxicity Studies
Toxicity prediction was carried out based on the validated and constructed models in Discovery studio software [28,29] as follows: (i) FDA rodent carcinogenicity, which computes the probability of a submitted chemical structure being a carcinogen; (ii) carcinogenic potency TD50, which computes the tumorigenic dose rate 50 (TD50) of a compound in a chronic exposure toxicity test for a rodent [30]; (iii) rat maximum tolerated dose (MTD) of the examined compounds [31,32]; (iv) rat chronic LOAEL, which predicts the rat chronic lowest observed adverse effect level (LOAEL) of a compound [33,34]; (v) ocular irritancy, which predicts the incidence and severity of ocular irritation of a particular compound in the Draize test (Wilhelmus, 2001 #41); and (vi) skin irritancy, which predicts the incidence and severity of the skin irritation of a particular compound in a rabbit skin irritancy test [35].
As shown in Table 4, the tested compounds showed in silico low toxicity against the tested models. For the FDA rodent carcinogenicity model, all examined compounds were expected to be non-carcinogenic.
For the carcinogenic potency TD50 rat model, the compounds showed TD50 values ranging from 4.059 to 5.463 mg/kg body weight/day. Such values are almost equal to the reference drug, sofosbuvir (4.929 mg/kg body weight/day). Regarding the rat maximum tolerated dose model, the compounds showed a maximum tolerated dose with a range of 0.021 to 0.063 g/kg body weight, which is almost equal to that of sofosbuvir (0.028 g/kg body weight). For the rat chronic LOAEL model, the tested compounds showed LOAEL values ranging from 0.0011 to 0.0026 g/kg body weight. These values are almost equal to or higher than sofosbuvir (0.0011 g/kg body weight). Moreover, all compounds were predicted to have a mild and moderate irritant effect against the skin irritancy and ocular irritancy models, respectively. Figure 13. ADMET_AlogP98: lipid-water partition coefficient; ADMET_PSA_2D: polar molecular surface area. Each drug is plotted with a 2D polar surface area (PSA_2D) against its computed atom-type partition coefficient (ALogP98). The area encompassed by the ellipse represents good absorption without any violation of the ADMET properties. The 95% and 99% confidence limit ellipses corresponding to the blood-brain barrier (BBB) and intestinal absorption models are indicated.

Toxicity Studies
Toxicity prediction was carried out based on the validated and constructed models in Discovery studio software [28,29] as follows: (i) FDA rodent carcinogenicity, which computes the probability of a submitted chemical structure being a carcinogen; (ii) carcinogenic potency TD 50, which computes the tumorigenic dose rate 50 (TD 50 ) of a compound in a chronic exposure toxicity test for a rodent [30]; (iii) rat maximum tolerated dose (MTD) of the examined compounds [31,32]; (iv) rat chronic LOAEL, which predicts the rat chronic lowest observed adverse effect level (LOAEL) of a compound [33,34]; (v) ocular irritancy, which predicts the incidence and severity of ocular irritation of a particular compound in the Draize test (Wilhelmus, 2001 #41); and (vi) skin irritancy, which predicts the incidence and severity of the skin irritation of a particular compound in a rabbit skin irritancy test [35].
As shown in Table 4, the tested compounds showed in silico low toxicity against the tested models. For the FDA rodent carcinogenicity model, all examined compounds were expected to be non-carcinogenic.
For the carcinogenic potency TD 50 rat model, the compounds showed TD 50 values ranging from 4.059 to 5.463 mg/kg body weight/day. Such values are almost equal to the reference drug, sofosbuvir (4.929 mg/kg body weight/day). Regarding the rat maximum tolerated dose model, the compounds showed a maximum tolerated dose with a range of 0.021 to 0.063 g/kg body weight, which is almost equal to that of sofosbuvir (0.028 g/kg body weight). For the rat chronic LOAEL model, the tested compounds showed LOAEL values ranging from 0.0011 to 0.0026 g/kg body weight. These values are almost equal to or higher than sofosbuvir (0.0011 g/kg body weight). Moreover, all compounds were predicted to have a mild and moderate irritant effect against the skin irritancy and ocular irritancy models, respectively.

DFT Studies
The calculated values of the global chemical reactivity parameters and the free energy and the dipole moment for the parent nucleosides of the active phosphoramidate prodrugs are shown in Table 5. Compound 4 was excluded from the DFT studies, which focused only on 2 -halogenated nucleosides. The chemical reactivity values vary with different halogen substituents and their stereochemistry. The electronic transition from the HOMO to LUMO represents the energy gap (E gap ), which describes the molecular reactivity. The value of the energy gap between HOMO and LUMO is directly proportional to the measure of molecules' excitability; the smaller Egap in a molecule, the easier its electronic excitement between HOMO and LUMO and vice versa [36]. Moreover, a lower E gap value would exhibit the eventual charge transfer interaction, which is generally responsible for the molecular bioactivity. Compounds 9 and 13 represent the highest Egap value (5.61 eV). Conversely, a larger E gap provides lower chemical reactivity. The LUMO energy describes the electron affinity (EA) as the amount of energy released when an electron is added to the molecule to form a negative ion. A higher EA of a molecule suggests an electron acceptor in a charge-transfer reaction while a molecule with a lower EA is more likely to be an electron donor. Thus, 3, with the lowest EA (1.48), has a higher electron donor capacity while 6, representing the highest EA (1.98), has a lower donor capacity. Inversely, the energy of the HOMO describes the ionization potential (IP). The electrophilicity (ω) is an index used to quantify the global electrophilic nature, meaning that 7, which has the lowest ω (3.48), is the most electrophile. This index also measures the stabilization in energy when the system acquires an additional electronic charge; the smaller the (ω) value, the stronger the molecule's stability.
The steric clash between 2 -substituted halogenated nucleoside analogs plays a pivotal role in the binding interaction of the 2 -β substitution and the active site in the target enzyme [18,19]. In order to understand the relationship between the van der Waals radii and volumes of 2 -substituted halogenated nucleosides and their biological activities, the Mulliken charges and bond length were calculated using DFT at the B3LYP/6-311+g(d,p) level of theory in addition to the Bondi radii and their volumes (Table 6) [37,38]. The results clearly indicate the bioisosteric effect between the methyl, chloro, and bromo groups, leading to similar biological activity. Table 6. Mulliken charges' distribution, van der Waals radii, and volumes of 2 -substituted nucleosides. Several studies were conducted in order to correlate different biological activities with chemical reactivity parameters such as the hardness, softness, chemical potential, electronegativity, electrophilicity, and other electronic parameters [36]. Simple linear regression of the biologically active and inactive parent nucleosides of the phosphoramidate prodrugs vs. global reactivity parameters of HCV genotypes 1a and 1b is shown in Table 7, representing the Pearson correlation coefficient (r) and the coefficient of determination (r 2 ). The significant statistical relationship (p < 0.05) of the Pearson coefficient between the global chemical reactivity parameters (descriptors) and the biological activity (predictors) against HCV genotype 1a using the B3LYP/6-311+g(d,p) level of theory was obtained for E gap , global electron affinity, chemical potential, electrophilicity, hardness, and softness with an insignificant p-value (>0.05) for HCV GT1a and GT1b. However, the results were slightly improved for HCV genotype 1a and 1b with the global ionization potential (r I = 0.67 and 0.72) and electronegativity (rχ = 0.63 and 0.63) with a significant p-value (<0.05) for HCV GT1a and GT1b. the Pearson coefficient was positive and directly proportional to the biological activity. Multiple linear regression for statistically significant global reactivity indexes of the ionization potential and electronegativity was achieved and improved the Pearson correlation coefficient: r = 0.72 and r = 0.75 for GT1a and GT1b, respectively. The parent 12 halogenated nucleosides were geometrically optimized at the groundstate using the density functional theory (DFT) at the B3LYP/6-311+g(d,p) level of theory, and MEP maps were calculated ( Figure 14). The molecular electrostatic potential (MEP) is the potential that a positively charged unit would experience at any point surrounding the molecule due to the electron density distribution. The electrostatic potential may predict the chemical reactivity of molecules since regions of negative potential are expected to be sites of protonation and nucleophilic attack while regions of positive potential may indicate electrophilic sites. to the biological activity. Multiple linear regression for statistically significant global reactivity indexes of the ionization potential and electronegativity was achieved and improved the Pearson correlation coefficient: r = 0.72 and r = 0.75 for GT1a and GT1b, respectively. The parent 12 halogenated nucleosides were geometrically optimized at the groundstate using the density functional theory (DFT) at the B3LYP/6-311+g(d,p) level of theory, and MEP maps were calculated ( Figure 14). The molecular electrostatic potential (MEP) is the potential that a positively charged unit would experience at any point surrounding the molecule due to the electron density distribution. The electrostatic potential may predict the chemical reactivity of molecules since regions of negative potential are expected to be sites of protonation and nucleophilic attack while regions of positive potential may indicate electrophilic sites.

Molecular Dynamics Simulation Studies
The molecular dynamics (MD) simulation experiment has an extra advantage over molecular docking, in addition to its ability to predict the binding mode correctly. It can predict the conformational changes of the ligand and the protein after the binding process [39]. Additionally, MD studies can accurately compute several parameters related to the change in the energy of the protein-ligand complex for a previously determined time. Accordingly, it precisely characterizes the stability, binding mode, and flexibility of the active compound inside the target enzyme [40]. The first MD simulation study of an enzyme was published in the Nature journal in 1977 [41]. Luckily, due to the recent supercomputer advancements, particularly the modern graphics processing parts, MD simulations analysis has become much more approachable, powerful, and precise [42]. The molecular docking and molecular dynamics simulations were performed for quinolinone derivatives with NS5B as a useful tool in order to design potent anti-HCV molecules [43].
Herein, the dynamics, together with the conformational variations of the HCV-NS5B GT1a-compound 6 complex, were calculated as a root mean square deviation (RMSD) over 100 ns. At first, the triphosphate groups were completely protonated, and Mg ++ was removed. The results indicate that the HCV-NS5B GT1a, compound 6, and HCV-NS5B GT1a-compound 6 complex displayed lower RMSD values, indicating greater stability. Although the HCV-NS5B GT1a-compound 6 complex showed some fluctuations, it revealed good stability after 90 ns~ ( Figure 15A). The flexibility of the HCV-NS5B GT1a residues was estimated in terms of the root mean square fluctuation (RMSF) to discover the fluctuated regions due to binding. As shown in Figure 15B, the binding of compound 6 did not cause major fluctuations in HCV-NS5B GT1a. Furthermore, the radius of gyration (Rg) was figured as a pivotal parameter associated with the protein constancy and stability [44,45]. The Rg of the HCV-NS5B GT1a-compound 6 complex ( Figure 15C) was computed over 100 ns to be lower at the end of the study than the starting period. Additionally, the interaction between HCV-NS5B GT1a-compound 6 complex and the encircling solvents was estimated by the solvent accessible surface area (SASA) over 100 ns. SASA of the HCV-NS5B GT1a-compound 6 complex was calculated to explore the conformational changes that occurred due to binding. Interestingly, HCV-NS5B GT1a demonstrated a decrease in the SASA value through the 100 ns of the study ( Figure 15D). The MD simulation studies revealed that the highest number of hydrogen bond conformations of the HCV-NS5B GT1a-compound 6 complex that was formed was up to 10 hydrogen bonds ( Figure 15E).

Molecular Dynamics Simulation Studies
The molecular dynamics (MD) simulation experiment has an extra advantage over molecular docking, in addition to its ability to predict the binding mode correctly. It can predict the conformational changes of the ligand and the protein after the binding process [39]. Additionally, MD studies can accurately compute several parameters related to the change in the energy of the protein-ligand complex for a previously determined time. Accordingly, it precisely characterizes the stability, binding mode, and flexibility of the active compound inside the target enzyme [40]. The first MD simulation study of an enzyme was published in the Nature journal in 1977 [41]. Luckily, due to the recent supercomputer advancements, particularly the modern graphics processing parts, MD simulations analysis has become much more approachable, powerful, and precise [42]. The molecular docking and molecular dynamics simulations were performed for quinolinone derivatives with NS5B as a useful tool in order to design potent anti-HCV molecules [43].
Herein, the dynamics, together with the conformational variations of the HCV-NS5B GT1a-compound 6 complex, were calculated as a root mean square deviation (RMSD) over 100 ns. At first, the triphosphate groups were completely protonated, and Mg ++ was removed. The results indicate that the HCV-NS5B GT1a, compound 6, and HCV-NS5B GT1a-compound 6 complex displayed lower RMSD values, indicating greater stability. Although the HCV-NS5B GT1a-compound 6 complex showed some fluctuations, it revealed good stability after 90 ns~( Figure 15A). The flexibility of the HCV-NS5B GT1a residues was estimated in terms of the root mean square fluctuation (RMSF) to discover the fluctuated regions due to binding. As shown in Figure 15B, the binding of compound 6 did not cause major fluctuations in HCV-NS5B GT1a. Furthermore, the radius of gyration (Rg) was figured as a pivotal parameter associated with the protein constancy and stability [44,45]. The Rg of the HCV-NS5B GT1a-compound 6 complex ( Figure 15C) was computed over 100 ns to be lower at the end of the study than the starting period. Additionally, the interaction between HCV-NS5B GT1a-compound 6 complex and the encircling solvents was estimated by the solvent accessible surface area (SASA) over 100 ns. SASA of the HCV-NS5B GT1acompound 6 complex was calculated to explore the conformational changes that occurred due to binding. Interestingly, HCV-NS5B GT1a demonstrated a decrease in the SASA value through the 100 ns of the study ( Figure 15D). The MD simulation studies revealed that the highest number of hydrogen bond conformations of the HCV-NS5B GT1a-compound 6 complex that was formed was up to 10 hydrogen bonds ( Figure 15E

MM-PBSA Studies
The MM-PBSA is an advanced, reliable tool used to estimate the free binding energy of a compound and protein. MM-PBSA is favored over other approaches such as the free energy perturbation and the thermodynamic integration because it is simpler, faster, and produces more consistent results [46]. The binding free energy of the HCV-NS5B GT1acompound 6 complex was estimated in the final 20 ns of the MD study at an interval of 100 ps from the MD trajectories. Compound 6 demonstrated a low binding free energy of -1663 KJ/mol with HCV-NS5B GT1a ( Figure 16A). Furthermore, the total binding energy of the HCV-NS5B GT1a-compound 6 complex was fragmented to discover the share of every amino acid residue of the HCV-NS5B GT1a in the binding with compound 6.

MM-PBSA Studies
The MM-PBSA is an advanced, reliable tool used to estimate the free binding energy of a compound and protein. MM-PBSA is favored over other approaches such as the free energy perturbation and the thermodynamic integration because it is simpler, faster, and produces more consistent results [46]. The binding free energy of the HCV-NS5B GT1acompound 6 complex was estimated in the final 20 ns of the MD study at an interval of 100 ps from the MD trajectories. Compound 6 demonstrated a low binding free energy of −1663 KJ/mol with HCV-NS5B GT1a ( Figure 16A). Furthermore, the total binding energy of the HCV-NS5B GT1a-compound 6 complex was fragmented to discover the share of every amino acid residue of the HCV-NS5B GT1a in the binding with compound 6.

Docking Studies
The crystal structures of the target proteins: (i) non-structural protein 5B (NS5BGT-1a) (PDB ID: 4khm, resolution: 1.70 Å) and (ii) non-structural protein 5B (NS5BGT-1b) (PDB ID: 4kai, resolution: 2.30 Å), were obtained from Protein Data Bank (http://www.pdb.org, accessed on 18 November 2020). Molecular Operating Environment (MOE) was utilized for the docking analysis [47]. The free binding energies and the binding modes of the tested compounds against the protein 5B were discussed. Firstly, the water molecules were removed from protein 5B's crystal structure, leaving one chain, which is vital for the binding. IPV, the co-crystallized ligand, was utilized as a reference ligand. Next, protein 5B's structure was protonated, and the hydrogen atoms were hidden. Additionally, all the triphosphate groups were completely protonated without any charge. Then, the binding pockets of the protein 5B were defined, and the energy was minimized [47,48].
The structures of the tested molecules and IPV ( Figure 10) were drawn using Chem-BioDraw Ultra 14.0 and kept as SDF files. Then, they were opened utilizing MOE software,

Docking Studies
The crystal structures of the target proteins: (i) non-structural protein 5B (NS5BGT-1a) (PDB ID: 4khm, resolution: 1.70 Å) and (ii) non-structural protein 5B (NS5BGT-1b) (PDB ID: 4kai, resolution: 2.30 Å), were obtained from Protein Data Bank (http://www.pdb.org, accessed on 18 November 2020). Molecular Operating Environment (MOE) was utilized for the docking analysis [47]. The free binding energies and the binding modes of the tested compounds against the protein 5B were discussed. Firstly, the water molecules were removed from protein 5B's crystal structure, leaving one chain, which is vital for the binding. IPV, the co-crystallized ligand, was utilized as a reference ligand. Next, protein 5B's structure was protonated, and the hydrogen atoms were hidden. Additionally, all the triphosphate groups were completely protonated without any charge. Then, the binding pockets of the protein 5B were defined, and the energy was minimized [47,48].
The structures of the tested molecules and IPV ( Figure 10) were drawn using Chem-BioDraw Ultra 14.0 and kept as SDF files. Then, they were opened utilizing MOE software, and the 3D structures were protonated. Next, the compound's energy was minimized.
Long-range electrostatic interactions were treated using the particle-mesh Ewald (PME) method, where a grid spacing of 1.0 Å was used for all simulation cells. All covalent bonds involving hydrogen atoms were constrained using the SHAKE algorithm for consistency, and the same protocol for all MD simulations was applied.

Conclusions
Twelve 2 -halogenated nucleotide analogs with a potential HCV polymerase inhibitory effect were subjected to many computational techniques, including docking, ADMET, toxicity, QSAR, DFT, and molecular dynamics simulation studies. Docking studies were performed against HCV GT1a and GT1b using sofosbuvir as a reference compound. The tested compounds exhibited a good binding mode against the target receptors, with a mild effect for the substitution of the ribose sugar on the binding potential. In addition, it was found that the high hydrophobic group can produce a higher binding affinity. The ADMET studies revealed that the tested compounds exhibited a low BBB penetration level, an optimal range of aqueous solubility, good intestinal absorption levels, a non-inhibitory effect on CYP2D6, and a plasma protein-binding ability of less than 90%. Moreover, the toxicity studies exhibited that the tested compounds were expected to be non-carcinogenic and equal to sofosbuvir regarding the carcinogenic potency TD 50 , maximum tolerated dose, and LOAEL values. Compound 6 was further subjected to molecular dynamics (MD) simulation to explain the protein-ligand complex stability and contacts. With these encouraging results, all the compounds can be further explored for structural modification and detailed investigations to arrive at possibly newer potent anti-HCV agents with better therapeutic activity.