key: cord-0761484-zsr5zz2k
authors: Hagar, Mohamed; Ahmed, Hoda A.; Aljohani, Ghadah; Alhaddad, Omaima A.
title: Investigation of Some Antiviral N-Heterocycles as COVID 19 Drug: Molecular Docking and DFT Calculations
date: 2020-05-30
journal: Int J Mol Sci
DOI: 10.3390/ijms21113922
sha: 4ce66b84c868384afe19e92dc095a950f1b17b12
doc_id: 761484
cord_uid: zsr5zz2k

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 M(pro), 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.

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; 

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 (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.

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.

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.

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. 

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 

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 

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.

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.

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] .

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. 

Serial interval of novel coronavirus (COVID-19) infections

Emerging understanding of etiology and epidemiology of the novel coronavirus (COVID-19) infection in Wuhan, China

China Novel Coronavirus I, Research T. A novel coronavirus from patients with pneumonia in China

A pneumonia outbreak associated with a new coronavirus of probable bat origin

COVID-19: A new challenge for human beings

In vitro antiviral activity and projection of optimized dosing design of hydroxychloroquine for the treatment of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2)

Hydroxychloroquine, a less toxic derivative of chloroquine, is effective in inhibiting SARS-CoV-2 infection in vitro

Physiologically based pharmacokinetic modeling in regulatory science: An update from the US Food and Drug Administration's Office of Clinical Pharmacology

New insights on the antiviral effects of chloroquine against coronavirus: What to expect for COVID-19?

Therapeutic efficacy of the small molecule GS-5734 against Ebola virus in rhesus monkeys

Remdesivir and chloroquine effectively inhibit the recently emerged novel coronavirus (2019-nCoV) in vitro

The $2.6 billion pill-methodologic and policy considerations

In silico oncology drug repositioning and polypharmacology

Individualized network-based drug repositioning infrastructure for precision oncology in the panomics era

Drug repurposing: New treatments for zika virus infection?

Network-based drug repurposing for novel coronavirus 2019-nCoV/SARS-CoV-2

Research and development on Therapeutic Agents and Vaccines for COVID-19 and Related Human Coronavirus Diseases

Therapeutic opportunities to manage COVID-19/SARS-CoV-2 infection: Present and future

Potential interventions for novel coronavirus in China: A systemic review

COVID-19: An update on the epidemiological, clinical, preventive and therapeutic evidence and guidelines of integrative Chinese-Western medicine for the management of 2019 novel coronavirus disease

Comparative therapeutic efficacy of remdesivir and combination lopinavir, ritonavir, and interferon beta against MERS-CoV

Unrevealing sequence and structural features of novel coronavirus using in silico approaches: The main protease as molecular target

Computational network biology: Data, model, and applications

In silico and in vitro analysis of small molecules and natural compounds targeting the 3CL protease of feline infectious peritonitis virus

N-alkyl 2-pyridone versus O-alkyl 2-pyridol: Ultrasonic synthesis, DFT, docking studies and their antimicrobial evaluation

A new chiral Boron-dipyrromethene (BODIPY)-based fluorescent probe: Molecular docking, DFT, antibacterial and antioxidant approaches

Computational and molecular docking approaches of a New axially chiral BODIPY fluorescent dye

Heterocyclic carbene metal complexes in medicinal chemistry

Heterocyclic naphthalimides as new skeleton structure of compounds with increasingly expanding relational medicinal applications

Favipiravir (T-705), a novel viral RNA polymerase inhibitor

Systematic review of amodiaquine treatment in uncomplicated malaria

The effect of 2 -fluoro-2 -deoxycytidine on herpes virus growth

The clinical pharmacology of ribavirin

Specific interactions between tau protein and curcumin derivatives: Molecular docking and ab initio molecular orbital simulations

DFT investigations of linear Zn3-type complex with compartmental N/O-donor Schiff base: Synthesis, characterizations, crystal structure, fluorescence and molecular docking

Spectroscopic identification, structural features, Hirshfeld surface analysis and molecular docking studies on stiripentol: An orphan antiepileptic drug

NMR and electronic spectral study on 1-phenyl-1-propanol by quantum computational methods

Modulating hydrogen-bond basicity within the context of protein-ligand binding: A case study with thrombin inhibitors that reveals a dominating role for desolvation

Role of resultant dipole moment in mechanical dissociation of biological complexes

Molecular recognition of azelaic acid and related molecules with DNA polymerase I investigated by molecular modeling calculations

Physiochemical, molecular docking, and pharmacokinetic studies of Naproxen and its modified derivatives based on DFT

Gaussian 09, Revision a. 02

AutoDock4 and AutoDockTools4: Automated docking with selective receptor 828 flexibility

This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license

Funding: This research received no external funding.

The authors declare no conflict of interest.