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

1Novel 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. Coronaviruses belong to the family of Coronaviridae, which is a family of enveloped- single 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. COVID-19 is a novel coronavirus of the β genus; it is round in shape or crown-shaped and has a diameter of ≈ 60–140 nm, as per observations under the electron microscopy. The incubation period of COVID-19 ranges from 2 to 14 days. According to the World Health Organization, there have been around 1 million and 550 thousand reported cases of COVID-19 disease and more than 90,000 reported deaths to date (11/04/2020). Moreover, SARS-COV-2, the big family of COVID-19, is sensitive to heat and ultraviolet (UV) rays. Furthermore, they are inactivated at 56 °C for half hour and can be stored at −80 °C for many years. Additionally, 75% ethanol, chlorine containing disinfectants and peracetic acid can effectively inactivate COVID-19. In general, the major principle to prevent and control the infectious disease is eliminating the source of infections thereby protecting the very sensitive population. The hydroxychloroquine an approved drug for malaria disease by FDA was explored as a medication for SARS-CoV-2. Previous reports revealed that, the chloroquine and hydroxychloroquine can inhibit the coronavirus (COVID-19) by changing the pH at the surface of the cell membrane. This action can inhibit the attachment of the virus to the cell membrane. In addition, it can prevent nucleic acid replication, glycosylation of viral proteins, virus assembly, new virus particle delivery, virus release, and other mechanisms to obtain its antiviral effects. Remdesivir, Ebola’s medication demonstrated activity against the coronavirus family in 2017 and 2020. To date, there are no SARS-CoV-2 vaccines available, however several domestic and foreign research institutions and enterprises have used several methods, including mRNA nano-vaccine technology, recombinant or inactivated vaccine technology, and DNA vaccine technology, to develop coronavirus vaccine. The treatment of this virus has not been achieved by any available antiviral agent. Hence, our present work suggests that selection of some derivatives with the appropriate viral restraining mechanisms can offer promising results.

2Although the development in the biomedical research has widespread area over the past two decades, the yearly number of new treatments agreed by the U.S. Food and Drug Administration (FDA) has remained relatively limited. Compared to de novo drugs and randomized clinical studies, the drug repurposing for virus including the novel coronavirus is represented as an impacted drug discovery strategy from existing drugs and could essentially shorten the time and reduce the cost. Recently, Chinese medicinal plants were classified as antiviral inhibitors of the novel COVID 19. The anti-coronavirus effects of these natural derivatives have been confirmed in vitro by direct loading onto cultured cells, however it does not guarantee their effectiveness in vivo. Some protease inhibitors may have an antiviral effect by blocking coronavirus main protease. This assumption was supported by FDA approval about the use of ritonavir/lopinavir, based on data obtained from in vitro studies. So, a recent comparative investigation was established to rationalize the potential use of protease inhibitors as a treatment against infections of SARS-CoV-2.

3Computational aspects offer new testable hypotheses for regular drugs involved for novel coronavirus. Recently, a virtual screening technique was carried out to identify the active site on the viral protease for the binding of many natural compounds through the molecular docking and cell-based assays.

4Recently, our research group has focused on the determination of the molecular geometry of synthesized materials via relating the desired properties from experimental evaluation with the estimated parameters from computational calculations.

5Heterocycles are widely investigated for possible medicinal applications. Favipiravir (1), amodiaquine (2), 2′-fluoro-2′-deoxycytidine (3), and ribavirin (4) are known as antiviral drugs. They inhibit various RNA and DNA viruses. Hence, we aim to determine whether the protease of COVID-19 can be a target protein of these nucleotides in silico using molecular docking. Moreover, a comparative study between these drugs with the FDA approved remdesivir and hydrocloroquine antiviral drugs against a broad range of RNA viruses has been established to investigate the effectiveness of the drugs as inhibitors for COVID 19.

6The 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.

7The 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 Bohr3), however, drug 2 of the highest complexity is predicted to have the highest polarizability, 268 Bohr3.

8Frontier 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.

9The 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.

10Recently, many reports showed that the FMOs have to be taken into consideration in investigation of the structure activity relationships. 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.

11The EHOMO and ELUMO are indicators for the prediction of the ionization potential (I= −EHOMO) and the electron affinity (A-ELUMO) 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:

12The χ 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. 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.

13The 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.

14To 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.

15The 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.

16The 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 drug 2 is recognized on C22 and C27 sites. However, N4 and N3 are the most negative charges of drugs 3 and 4 and their respective positively charged atoms are C5 and C15. The positively charged centers are the most susceptible sites for nucleophilic attacks i.e., electron donation. However, the most negatively charged centers are the most susceptible sites for electrophilic one.

17Docking simulation studies were done to predict the binding mode of drugs 1–4 at the SARS-CoV-2 Mpro 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 Mpro 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 (N11, N31, and O29) 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.

18As 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 N11 (−0.729) N31 (−0.742) rather than O29 (−0.612) for drug 2. However, the negative electrostatic potential is localized on N1 (−0.336) N3 (−0.713) rather than O7 (−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.

19On 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. Another factor that could affect the degree of the interaction of these drugs with the protein is the dipole moment. 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.

20The theoretical estimations were carried out for the investigated drugs by Gaussian 09 software, (2009, Gaussian 09, revision a. 02, gaussian. Inc.: Wallingford, CT, USA). 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. 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.

21Docking studies were carried out using the AUTODOCK 4.2 program. 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).

22Before 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.

23Four 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 N11, N31, and O29 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.