Molecular Docking and Conceptual DFT-Based Study of Some Potential SARS-CoV-2 Inhibitors

Nowadays, the main effort of the scientific community is focused on the search of specific drugs for the inhibition of the Severe Acute Respiratory Syndrome—Coronavirus 2 (SARS-CoV-2), which is responsible for the Coronavirus Disease 19 or COVID-19. With the same objective in mind, a Molecular Docking study was performed in this work in order to discover information about some antiviral drugs of common use as protease inhibitors. As a complement of this research, a chemical reactivity study of these potential drugs was done with the aim of finding a relationship between the inhibition ability and the chemical reactivity. The results presented in this research con-stitute one of the first predictions aimed to identify the best potential compounds for this purpose while at the same time verifying the validity of the employed theoretical and computational methodology. By means of the analysis of the number of hydrogen bonds as well as the binding energies coming from the Molecular Docking study, it can be said that Telaprevir, Nelfinavir and Indinavir have the highest probability of success as potential inhibitors of SARS-CoV-2.


Introduction
During the rapid development of the COVID-19 disease, researchers over the whole world are hard-working in order to develop a successful therapy. The main protease proteins of SARS-CoV-2 have been identified as the drug targets required to process the translation of the proteins from the RNA virus. In particular, the main protease (Mpro, also known as 3CLpro) is one of the coronavi-DOI: 10 [3]. Thus, inhibition of Mpro would prevent the virus from replication and therefore constitutes one of the potential anticoronaviral strategies [1] [2] [3] [4]. For these reasons, the main objective of some of recent investigations has been finding adequate drugs with inhibitory activity properties of this particular protease protein [3]- [9].
Viral proteases show differences depending on the type of virus that express them. Among other aspects, they may vary in their substrate specificity as well as in their amino acid sequences and their subunits. However, all of them have a series of regions of substrate specificity where they can act as receptors, the so-called active center [10]. In order to find the particular region of one of the main proteases of SARS CoV-2 as well as the interactions with the amino acids present in these areas, a Molecular Docking process with molecules with potential inhibition activity can be developed with the aim of finding information about the protein-substrate interactions. With the objective of identifying if already well-known drugs can be used to treat this new disease that nowadays is affecting more than 200 countries and territories around the world [11], a series of protease inhibitors that could block the activity of the SARS-CoV-2 protease enzyme can be tested by using Molecular Docking strategies. These kinds of systems work by blocking the division of the large polyproteins into the smaller fragments that are necessary for the assembling of new viral particles [1] [2] [3] [4].
In the drug discovery research, the identification and validation of lead compounds and the determination of active binding sites of biological targets related to a certain lead compound performed through wet lab experiments are quite expensive and time consuming [12]. Computational determinations reduce effectively the time required to obtain useful drugs and decrease their associated economic costs, making it possible the proposal of new potential drugs with low expenditures and selective targeting [13].
Considering these facts, the overall goal of this research is the discovery of significant binding affinities between some protease inhibitors and a selected protein chosen as the receptor, considering it in a rigid conformation and docking the protease inhibitors into its active site. Once the enzyme-substrate interaction is selected, the binding site at several discrete points is tested and the interaction energies between the molecules and the atoms within the receptor are calculated. This information is used to identify energetically favourable and unfavourable regions for specific ligand-receptor interactions that could guide a correct pharmacophore feature placement or aid in the ligand design and optimization [14] [15] [16].
The selection of the protease inhibitors to be tested in this research was inspired by the recently experimental evaluation of their antiviral activity by Yamamoto and coworkers [17].  [27], while the crystallographic structure of the SARS-CoV-2 main protease was obtained from Protein Data Bank (PDB) with code 6M03 [28]. The physical and chemical properties of these molecules as well as their pharmacokinetics are readily available through the free online PubChem database (https://pubchem.ncbi.nlm.nih.gov/). A graphical representation of the two-dimensional structures of the considered molecules is displayed in Figure 1.
In a complementary way, the chemical reactivity properties of the potential inhibitors were determined through a calculation of the reactivity descriptors that arise from Conceptual Density Functional Theory (CDFT) [29] [30] [31] [32]. The chemical reactivity descriptors affect the selectivity of the molecular system and play an important role in the prediction of the electrodonating or electroaccepting character of the systems [32]- [39]. The electronic configuration details help to describe the recognition process under negative and positive charges, exposing those interactions that can be difficult to identify with non bonded interactions [40]. The use of the concepts derived from this theory as well as their associated descriptors allowed us to find the reactivity of the molecular systems in terms of their electronic densities which was an additional contribution in the validation of the lead compounds. The determination of the Conceptual DFT-based chemical reactivity descriptors of the potential SARS-CoV-2

Molecular Systems Relaxation and Optimization
The molecular structures of the ligands were obtained from the online ZINC15 database [48]. Their structures were geometrically optimized and the most stable conformer in each case was selected. Random sampling was used in the process After the conformational determination procedure, the selected molecular structures were geometrically reoptimized followed by the analysis of the vibrational frequencies using the Gaussian 09 series of programs [51]. The procedure was performed in the context of the Density Functional Theory (DFT) methodology using the MN12SX density functional [52]

KID Procedure
It is usually assumed that the goodness of a given density functional can be esti- + . An additional descriptor ΔSL has been designed [41] [42] [43] to help in the verification of the accuracy of the KID approximation by comparing the HOMO energy of the radical anion with the energy of the LUMO of the neutral species. Although the Koopmans-complaint behavior of the MN12SX density functional has been proven previously for the case of peptides [41] [42] [43], we think that it is worth to perform a further validation for the case of the molecules considered in the present study.

Chemical Reactivity Properties
On the basis of the molecular structures obtained from the geometry optimizations, the calculation of the KID parameters as well as the Conceptual DFT-based chemical reactivity descriptors was pursued by resorting to our in-house developed "CDFT" application [41] [42] [43]. The parameters found were the ionization potential (I) and the electron affinity (A), directly from the results of the ground state computed energy, by using the same density functional and basis set as for the optimization step. This was followed by the determination of the global reactivity descriptors, which are: global hardness (η ), electronegativity ( χ , and electrophilicity ( ω ) [29] [30] [31]. Also, the electrodonating power ( ω − ) and electroaccepting power ( ω + ) [39] and finally, the net electrophilicity ( ω ± ∆ ) [56], whose meanings are described as follows: • Chemical Potential. This property measures the tendency of the electron to escape from systems in equilibrium. It is equal to minus the electronegativity. • Global Hardness. It quantifies the resistance of a molecule to intramolecular charge transfer.
• Electrophilicity Index. It is defined in terms of the square of chemical potential divided by the chemical hardness and it denotes the stabilization energy after a system gets an additional electronic charge from the external environment.
• Electrodonating Power. It represents the capability of a chemical system to donate a small fractional charge. • Electroaccepting Power. It represents the capability of a chemical system to accept a small fraction of charge.
• Net Electrophilicity. It is a measure of the relative electrophilicity.
Their mathematical definitions are expressed through the following equations [29] being H  and L  the HOMO and LUMO energies associated with each of the molecules considered in this work.

Molecular Docking
The crystal structure of COVID-19 main protease related to SARS-CoV-2 with PDB ID: 6M03 has been retrieved from the Protein Data Base (https://www.rcsb.org/structure/6M03) and was considered as the receptor in this study. The molecular docking was performed with the aid of AutoDock 4.2.6 [57]. As a part of the procedure, water molecules were eliminated and only polar hydrogens were added to the protein structure for the simulation process of the intermolecular interactions between the protein and the ligand in the active site.
The Lamarckian Genetic Algorithm (LGA) was considered for the estimation of the free energy change upon binding [58].

Geometry Optimization
The

KID Procedure Validation
Although the Koopmans-complaint behavior of the MN12SX density functional has been proven previously for the case of peptides [41] [42] [43], we think that it is worth to perform a further validation for the case of the molecules considered in the present study. This determination has been achieved by making use of the in-house developed CDFT software tool and the results of this analysis are shown in Table 1.
It can be seen from the results in Table 1 that the descriptors considered for the estimation of the goodness of the selected density functional through the KID procedure are very close to zero for all the studied molecules, providing an accurate justification for the choice of the MN12SX/Def2DZVP/H2O model chemistry employed for the computational determinations in this study.

Chemical Reactivity of the SARS-CoV-2 Inhibitors
With the most stable and relaxed molecular structure of the analyzed SARS-CoV-2 DOI: 10.4236/cmb.2020.104008 118 Computational Molecular Bioscience inhibitors, the analysis of their reactivity was performed using the descriptors mentioned above and included in Table 2.
The electronic affinities of inhibitors vary from 0.89 to 2.66 eV. The greater facility to form an anion is for Fipranavir. For the case of the ionization potential, Telaprevir is the system with the greatest potential of losing an electron with  [63] have proposed a Nucleophilicity index N through the consideration of the HOMO energy obtained through the KS scheme with an arbitrary shift of the origin taking the molecule of tetracyanoethylene (TCE) as a reference. On the basis of the previous definition and the scale established by these authors [60], it can be concluded that Nelfinavir can be regarded as a strong nucleophile because the value for the Nucleophilicity N is slightly greater than 3 eV, while the other inhibitors can be considered as moderate nucleophiles with the exception of Telaprevir with a N value of 1.77 eV, thus labeled as a marginal nucleophile.

Molecular Docking and Active Site Descriptors for the SARS-CoV-2 Inhibitors
The Molecular Docking calculations were performed for the interaction between the protease protein with PDB code 6M03 and twelve potential inhibitors of protease included in the ZINC15 database of the Food and Drug Administration (FDA). One of the important results that arise from the Molecular Docking process is the binding energy (BE) of the inhibitors. Negative values for this BE are an indication of the stability of the system and for the interaction with the active site of the protein. The docking procedure was evaluated for 10 poses.
The other important result from the molecular docking is the description of The active site of the macromolecule remains in the core of the protein structure for all the tested inhibitors. Some of the residues are linked together in protein sequences. The residues of the active site for each ligand are presented in Table 3.
There are residues repeated in the contour of the pose with the best score of each ligand. Lysine5 (LYS5) is present in the active site of eight residues: Amprenavir, Darunavir, Fosamprenavir, Indinavir, Nelfinavir, Saquinavir, Telaprevir

Hydrogen Bonds
The electrostatic interactions analysis between the ligands and the residues of the active site was performed with the aim of classifying them as strong, moderate or weak following the recommendation of Jeffrey [66]. It could be appreciated that seven of the twelve studied ligands presented hydrogen bonds which can be considered significative because the electrostatic interactions as hydrogen bonds are important as a measure of the highest probability of inhibitory effectiveness. This information is also summarized in Table 3.
• Amprenavir-Weak hydrogen bond between O7 of the (C=O) accepting group of the ligand and amino group of LYS5 and with bond distance of 2.24 Å. Moderate hydrogen bonds between carboxylic group of ASP289 and N10 from the amino of the ligand with bond distance of 2.10 Å.

Conclusions
In this research, we have performed the calculation of the chemical reactivity properties and a Molecular Docking study of the main protease proteins of SARS-CoV-2 6M03 with twelve known molecules with recognized inhibitory activity included in the ZINC15 drugs database of the Food and Drug Administration (FDA). Binding energies, hydrogen bonds and residues that form the active site were defined for all the ligands.
The results obtained from the determination of the chemical reactivity parameters indicate that the lowest values for the chemical hardness were for Tipranavir with 3.49 eV, Asunaprevir with 3.96 eV and Nelfinavir with 4.03 eV. This means that these ligands have the lowest resistance to change the shape of their electronic densities, and thus are the most reactive ones. According to the Binding Energy results obtained from the Molecular Docking procedure, Indinavir had the greatest stability with the studied protease followed by Lopinavir, Tipranavir and Nelfinavir with a small binding energy difference. The binding site within the biological target maintained the presence of some residues in simple form and linked in protein arrangements. LYS5, GLU290 and dimer LYS137-