Investigation of Some Antiviral N-Heterocycles as COVID 19 Drug: Molecular Docking and DFT Calculations

The novel coronavirus, COVID-19, caused by SARS-CoV-2, is a global health pandemic that started in December 2019. The effective drug target among coronaviruses is the main protease Mpro, because of its essential role in processing the polyproteins that are translated from the viral RNA. In this study, the bioactivity of some selected heterocyclic drugs named Favipiravir (1), Amodiaquine (2), 2′-Fluoro-2′-deoxycytidine (3), and Ribavirin (4) was evaluated as inhibitors and nucleotide analogues for COVID-19 using computational modeling strategies. The density functional theory (DFT) calculations were performed to estimate the thermal parameters, dipole moment, polarizability, and molecular electrostatic potential of the present drugs; additionally, Mulliken atomic charges of the drugs as well as the chemical reactivity descriptors were investigated. The nominated drugs were docked on SARS-CoV-2 main protease (PDB: 6LU7) to evaluate the binding affinity of these drugs. Besides, the computations data of DFT the docking simulation studies was predicted that the Amodiaquine (2) has the least binding energy (−7.77 Kcal/mol) and might serve as a good inhibitor to SARS-CoV-2 comparable with the approved medicines, hydroxychloroquine, and remdesivir which have binding affinity −6.06 and −4.96 Kcal/mol, respectively. The high binding affinity of 2 was attributed to the presence of three hydrogen bonds along with different hydrophobic interactions between the drug and the critical amino acids residues of the receptor. Finally, the estimated molecular electrostatic potential results by DFT were used to illustrate the molecular docking findings. The DFT calculations showed that drug 2 has the highest of lying HOMO, electrophilicity index, basicity, and dipole moment. All these parameters could share with different extent to significantly affect the binding affinity of these drugs with the active protein sites.


Introduction
Novel coronavirus (COVID-19) emerged as an infection and quickly spread to all countries and is primarily transmitted by contact with infected droplets of saliva or discharge from the nose when infected people cough or sneeze. The first identification of human coronaviruses was observed in the mid-1960s [1,2]. Coronaviruses belong to the family of Coronaviridae, which is a family of envelopedsingle stranded-positive sense RNA virus. In addition to, the family of Coronaviridae was divided into four genera: α, β, γ, and δ. Coronaviruses of α and β genera generally infect mammals and humans while the type of γ and δ genera mainly infect birds. That specification is according to the phylogenetic analysis and genome structure of coronaviruses [3]. COVID-19 is a novel coronavirus of the β genus;

DFT Calculations Studies
The theoretical DFT calculations were performed in gas phase by DFT method at B3LYP 6-311G (d,p) basis set. All optimum compounds (1-4, Figure 1) are stable, and this is approved in terms of the absence of the imaginary frequency. The results of the theoretical DFT calculations for all investigated drugs revealed the non-planarity except for compound 1. The estimated DFT calculations for thermal parameters, dipole moment and the polarizability of the drug derivatives 1-4 are summarized in Table 1 The theoretical DFT calculations were performed in gas phase by DFT method at B3LYP 6-311G (d,p) basis set. All optimum compounds (1-4, Figure 1) are stable, and this is approved in terms of the absence of the imaginary frequency. The results of the theoretical DFT calculations for all investigated drugs revealed the non-planarity except for compound 1. The estimated DFT calculations for thermal parameters, dipole moment and the polarizability of the drug derivatives 1-4 are summarized in Table 1. The DFT estimated data revealed that the dipole moment of the drugs under investigation is in the order of 2 ˃ 4 ˃ 1 ˃ 3. The high dipole moment 2 and 4 could illustrate their binding pose within a specific target protein and their results of the predicted binding affinity that will be discussed in the following molecular docking part. The polarizability of the materials depends on how the susceptibility of molecular system electron cloud be affected by approaching of a charge. Moreover, it depends on the complexity of the compounds as well as the size of the molecular structure. Molecules of the large size are more polarizable compounds. It is worth noting that the compound 1 is the smallest in size and has the least polarizability (76 Bohr 3 ), however, drug 2 of the highest complexity is predicted to have the highest polarizability, 268 Bohr 3 .   The DFT estimated data revealed that the dipole moment of the drugs under investigation is in the order of 2 <4 <1 <3. The high dipole moment 2 and 4 could illustrate their binding pose within a specific target protein and their results of the predicted binding affinity that will be discussed in the following molecular docking part. The polarizability of the materials depends on how the susceptibility of molecular system electron cloud be affected by approaching of a charge. Moreover, it depends on the complexity of the compounds as well as the size of the molecular structure. Molecules of the large size are more polarizable compounds. It is worth noting that the compound 1 is the smallest in size and has the least polarizability (76 Bohr 3 ), however, drug 2 of the highest complexity is predicted to have the highest polarizability, 268 Bohr 3 .

Frontier Molecular Orbitals
Frontier molecular orbitals (FMOs) are the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). The HOMO is the highest energy orbital occupied with electrons, so it is an electron donor, while, LUMO is the lowest energy orbital that has a space to accept electrons, so it is an electron acceptor. These orbitals control the mode of the interaction of the drugs with other molecules such as the interactions between these drugs and their receptors. The frontier molecular orbitals (FMO) can give realistic qualitative information about susceptibility of the electrons of the HOMO to transfer to the LUMO. Moreover, HOMO and LUMO are very important quantum chemical parameters to determine the reactivity of the molecules and are used to calculate many important parameters such as the chemical reactivity descriptors. The energies of the HOMOs and LUMOs of the studied compounds were calculated using DFT method at B3LYP 6-311G (d,p) basis set and are tabulated in Table 2. The isodensity surface plots of HOMO and LUMO for investigated compounds are shown in Figure 2.  The results of the FMOs energy analysis revealed that the energies of HOMOs of drugs 1 and 4 are lower compared with the other compounds 2 and 3. However, the destabilization of the LUMO level is found to be higher in 1 than the others. Consequently, the energy gap of studied drugs is in the order of 1 < 2< 3 < 4.
Recently, many reports showed that the FMOs have to be taken into consideration in investigation of the structure activity relationships [34][35][36]. The FMOs theory showed that the The results of the FMOs energy analysis revealed that the energies of HOMOs of drugs 1 and 4 are lower compared with the other compounds 2 and 3. However, the destabilization of the LUMO level is found to be higher in 1 than the others. Consequently, the energy gap of studied drugs is in the order of 1 < 2< 3 < 4.
Recently, many reports showed that the FMOs have to be taken into consideration in investigation of the structure activity relationships [34][35][36]. The FMOs theory showed that the energy level of the HOMO and the LUMO are the most significant aspects that impact the bioactivities of small structural drugs. Mainly HOMOs that offer electrons, however, the LUMOs accept electrons. Obviously, the level of energy of HOMOs are different for all investigated drugs. Compound 2 showed the most lying HOMO than the other drugs and consequently it could be a better electron donor drug. Interestingly, drug 4 of the largest energy gap ∆E = 5.54 eV, there are several hydrophilic interactions that could facilitate the binding with the receptors. This suggests that such hydrophilic interactions considerably impact the binding affinity of such small drugs to the receptors. The HOMO of a certain drug and the LUMO with the adjacent residues could share the orbital interactions during the binding process.

Chemical Reactivity Descriptors
The E HOMO and E LUMO are indicators for the prediction of the ionization potential (I= −E HOMO ) and the electron affinity (A-E LUMO ) of molecules. Besides the frontier molecular orbitals are used in estimation of other chemical reactivity descriptors such as electronegativity (χ), global hardness (η), softness (δ), and electrophilicity (ω). These are calculated according to the following equations [37]: The χ value is a prediction of the power of the molecule to attract electrons i.e., Lewis acid, while small values of (χ) are indication of a good base. The global hardness (η) is a degree of their charge transfer prohibition, however, the global softness (δ) characterizes the ability of a molecule to accept electrons [37]. Soft molecules are of a small energy gap between frontier molecular orbitals and are more reactive than the harder because they could easily transfer electrons to the acceptors.
The electrophilicity (ω), calculated from the electronegativity and chemical hardness, is an indicator of lower energy difference due to the highest electron movement between the acceptor, LUMO, and the donor, HOMO.

Molecular Electrostatic Potential (MEP)
To validate the evidence about the reactivity of the drug as inhibitors, the molecular electrostatic potential (MEP) is important to be calculated. Although the MEP gives an indication about the molecular size and shape of the positive, negative as well as the neutral electrostatic potential. These could be a tool to predict physicochemical property relationships with the molecular structure of the drugs under investigation. Moreover, the molecular electrostatic potential is a useful tool to estimate the reactivity of the drugs toward electrophilic and nucleophilic attacks.
The molecular electrostatic potential of the studied drugs (1-4) is calculated by the same method under the same base sets and seen in Figure 3. In the MEP, the maximum negative region is the preferred sites for electrophilic attack, indicated as red color. So, an attacking electrophile will be attracted by the negatively charged sites, and the opposite satiation for the blue regions. It is obvious that the molecular size and the shape as well as the orientation of the negative, positive, and the neutral electrostatic potential varied according to the drug because of the type of the atoms and its electronic nature. The difference in the mapping of the electrostatic potential around the drug could be principally responsible for variation of its binding affinity with the active sites receptor.

Molecular Electrostatic Potential (MEP)
To validate the evidence about the reactivity of the drug as inhibitors, the molecular electrostatic potential (MEP) is important to be calculated. Although the MEP gives an indication about the molecular size and shape of the positive, negative as well as the neutral electrostatic potential. These could be a tool to predict physicochemical property relationships with the molecular structure of the drugs under investigation. Moreover, the molecular electrostatic potential is a useful tool to estimate the reactivity of the drugs toward electrophilic and nucleophilic attacks.
The molecular electrostatic potential of the studied drugs (1-4) is calculated by the same method under the same base sets and seen in Figure 3. In the MEP, the maximum negative region is the preferred sites for electrophilic attack, indicated as red color. So, an attacking electrophile will be attracted by the negatively charged sites, and the opposite satiation for the blue regions. It is obvious that the molecular size and the shape as well as the orientation of the negative, positive, and the neutral electrostatic potential varied according to the drug because of the type of the atoms and its electronic nature. The difference in the mapping of the electrostatic potential around the drug could be principally responsible for variation of its binding affinity with the active sites receptor.

Mulliken Atomic Charges
The Mulliken atomic charges of the estimated drugs (1-4) were calculated the DFT using B3LYP 6 as a method at -311G (d,p) at a basis set, the data were tabulated in Table 3. It showed that the C7 is the most positive and O10 have the most negative charge for drug 1. On the other hand, it is observed that the most nucleophilic centers of drug 2 are N11 and N31 which are the most electrophilic susceptibility positions. On the other hand, it is obvious that the nucleophilic susceptibility of the

Mulliken Atomic Charges
The Mulliken atomic charges of the estimated drugs (1-4) were calculated the DFT using B3LYP 6 as a method at -311G (d,p) at a basis set, the data were tabulated in Table 3. It showed that the C 7 is the most positive and O 10

Molecular Docking
Docking simulation studies were done to predict the binding mode of drugs 1-4 at the SARS-CoV-2 M pro pocket. Similarly, the hydroxychloroquine and remdesivir were docked in the protein pocket as a control. The binding affinity of compounds 1-4 ranged between −4.06 to −7.77 Kcal/mol as presented in Table 4. However, the superpositions of the resulted poses of ligands 1-4 with the docking models of hydroxychloroquine and remdesivir, clearly showed that ligand 1 (violet stick) is not located in the catalytic dyad of the receptor (Figure 4). While, ligands 2-4 are located in the active pocket similar to hydroxychloroquine and remdesivir. Moreover, ligand 3 (green stick) did not demonstrate any significant interaction with the Cys-His catalytic dyad although it was oriented toward His-41 ( Figure 5). Strikingly, ligand 4 is predicted to bind to the catalytic dyad (Cys-145 and His-41) with strong hydrogen bonds along with other residues with binding energy (−4.69 Kcal/mol). The 2D representation of the binding mode of the ligand 4 in the receptor was predicted a high efficiency of the inhibitor with respect to the remdesivir as illustrated in Figure 5. Interestingly, ligand 2 with the least binding energy (−7.77 Kcal/mol) might be better inhibitor to SARS-CoV-2 M pro comparable with the hydroxychloroquine and remdesivir with binding affinity −6.06 and −4.96 Kcal/mol, respectively. Ligand 2 is interacted with the receptor with different hydrophobic interactions besides three hydrogen bonds with amino acids Leu141, Cys145, and Gln189 in a distance 2.07, 2.91, and 2.39 Å ( Figure 6). The molecular docking results revealed that the effective interactions of proteins with the drug 2 were found on the atoms (N 11 , N 31 , and O 29 ) and this could be attributed to the presence of lone pair electrons on these atoms. Moreover, the π-π stacking could be another hydrophobic interaction between the drug and the receptors. As previously discussed from the DFT calculations, the electrostatic potential gives an idea about the molecular size and shape of the positive, negative, as well as the neutral electrostatic potential to estimate the reactivity of the studied drugs toward electrophilic and nucleophilic attacks, moreover, van der Waals surface will provide some more negative electrostatic force found on N 11 (−0.729) N 31 (−0.742) rather than O 29 (−0.612) for drug 2. However, the negative electrostatic potential is localized on N 1 (−0.336) N 3 (−0.713) rather than O 7 (−0.490) for drug 4. These negative MEP charges could be used for H-bond formation or van der Walls interactions with protein receptors. It is obvious, that the negative electrostatic potential of drug 2 is higher than that of the drug 4, and this could be an illustration of the least binding affinity of 2 compared with other ligands 1, 3, and 4. effect on the binding affinity [11,38]. Another factor that could affect the degree of the interaction of these drugs with the protein is the dipole moment [39][40][41]. The calculated dipole moments of the investigated drugs are in the range of 14. 44-2.26 Debye in the order of 2 > 4 > 1> 3, Table 1. All these factors could share together with different extent to significantly impact the degree of the binding affinity of these drugs with the active protein sites. these drugs with the protein is the dipole moment [39][40][41]. The calculated dipole moments of the investigated drugs are in the range of 14.44-2.26 Debye in the order of 2 > 4 > 1> 3, Table 1. All these factors could share together with different extent to significantly impact the degree of the binding affinity of these drugs with the active protein sites. Underline: The most interactive amino acid in the binding pocket.   On the other hand, DFT calculations of the frontier molecular orbitals to make a comparison between the energy level of HOMO and LUMO as well as with their energy gap of the investigated drugs showed their impact on the bioactivity of these compounds. The energy levels of HOMOs are between −5.65 eV to −7.61 eV, however, the LUMO are in between −1.05 to 5.47 eV depending on the conjugation as well as the nature of substituents around the nucleus. Further, the order of the FMOs energy gap was 2.14, 4.00, 5.28, and 5.54 for 1, 2, 3, and 4, respectively. It is clear that, the drug 2 of the high lying HOMO is the most susceptible to be electron donor. Moreover, drugs 2 and 4 are considered good electrophiles because of their high electrophilicity indices 13.69 and 12.80, respectively. Moreover, drug 2 of a high basicity (χ = 3.70) rather than the others could be another effect on the binding affinity [11,38]. Another factor that could affect the degree of the interaction of these drugs with the protein is the dipole moment [39][40][41]. The calculated dipole moments of the investigated drugs are in the range of 14.44-2.26 Debye in the order of 2 > 4 > 1 > 3, Table 1. All these factors could share together with different extent to significantly impact the degree of the binding affinity of these drugs with the active protein sites.

DFT Calculations
The theoretical estimations were carried out for the investigated drugs by Gaussian 09 software, (2009, Gaussian 09, revision a. 02, gaussian. Inc.: Wallingford, CT, USA) [42]. B3LYP 6-311G basis set was nominated for the DFT calculations. The structural geometry was optimized by minimizing its energies compared to all geometrical variables without forcing any molecular symmetry restrictions. The molecular structure of the optimized drugs has been drawn by Gauss View [43]. The calculated frequency showed that all molecular structures of the investigated drugs was stationary points in the geometry optimization method without the presence of an imaginary frequency.

Molecular Docking Procedure
Docking studies were carried out using the AUTODOCK 4.2 program [44]. The crystal structure of COVID-19 main protease at 2.16 Å resolution was retrieved from the Protein Data Bank (PDB: 6LU7): https://www.rcsb.org/structure/6LU7. The 3D structures in PDB format of hydroxychloroquine (DB01611) and remdesivir (DB14761) were obtained from the DrugBank data base. The selected drugs 1-4 were extracted from the PubChem database in the SDF format and were converted to PDB format using PyMOL (The PyMOL Molecular Graphics System, v1.6-alpha; Schrodinger LLC, New York, NY, USA 2013).
Before docking compounds 1-4 on the target, the protein was edited using AutoDockTools (ADT). The water molecules were removed, the polar hydrogen atoms were added to the amino acid residues and Gasteiger charges were assigned to all atoms of the protein. Then the protein in PDBQT format was used as an input for the AUTOGRID program. AUTOGRID performed a pre-calculated atomic affinity grid maps for each atom type in the ligand plus an electrostatics map and a separate desolvation map present in the substrate molecule taking the entire protein as the search space. Flexible ligand docking was performed for each compound. Docking calculations were carried out using the Lamarckian genetic algorithm (LGA), and all parameters were the same for each docking. This process was repeated 100 times for each ligand, and the final mean affinity score was taken. The results were shown using Discovery Studio Visualizer v17.2.0.16349 [44].

Conclusions
Four known antiviral drugs, favipiravir (1), amodiaquine (2), 2 -fluoro-2 -deoxycytidine (3), and ribavirin (4), have been investigated as inhibitors for COVID-19 by DFT and molecular docking calculations. The results of the docking studies aimed at studying the binding mode of these drugs to the SARS-CoV 3CLpro. The results revealed that amodiaquine (2) and ribavirin (4) showed the best affinity with the target receptor, even though this behavior was not experimentally verified. It was predicted that the amodiaquine showed the lowest binding energy (−7.77 Kcal/mol) with respect to the other drugs. Moreover, amodiaquine binding affinity is lower than that of the approved medicines, hydroxychloroquine and remdesivir which have binding affinity −6.06 and −4.96 Kcal/mol; respectively. The results of the molecular docking have been illustrated in terms of the DFT calculations. The DFT results showed that amodiaquine (2) is the most lying HOMO, and consequently it could be the best to act as an electron donor. Moreover, the most electrophilic centers of the drug 2 are N 11 , N 31 , and O 29 and this could be attributed to the presence of the lone pair of electrons on these atoms. Further, the π-π stacking could be another hydrophobic interaction between the drug and the receptors. Additionally, a high basicity (χ = 3.70) and dipole moment (µ = 14.4 Debye) of drug 2 rather than the others, could be the other factors that enhanced the extent of the binding affinity. It could be concluded that these parameters share together with different magnitudes and affect the degree of the binding affinity of these drugs with the active protein sites to afford a certain degree of inhibition. Finally, from these in silico studies for drug 2, it is very promising to perform further in vitro and small animal model in vivo studies to establish a solid experimental evidence of its activity as inhibitors for COVID-19.