New findings on ligand series used as SARS-CoV-2 virus inhibitors within the frameworks of molecular docking, molecular quantum similarity and chemical reactivity indices [version 1; peer review: awaiting peer review]

Background: The severe acute respiratory syndrome coronavirus (SARS-CoV)-2 virus causes an infectious illness named coronavirus disease 2019 (COVID-19). SARS-CoV is a positive-sense single-stranded RNA virus from the Betacoronavirus genus. The SARS-CoV-2 RNA-dependent RNA polymerase (RdRp) has an important role in the viral life cycle and its active site is a very accessible region, thus a potential therapeutic approach may be to target this region to study the inhibition of viral replication. Various preexisting drugs have been proposed for the treatment of COVID-19 and the use of existing antiviral agents may reduce the time and cost of new drug discoveries, but the efficacy of these drugs is limited. Therefore, the aim of the present study was to evaluate a number of ligands used as SARS-CoV-2 virus inhibitors to determine the suitability of them for potential COVID-19 treatment. Methods: In this study, we selected a series of ligands used as SARS-CoV-2 virus inhibitors such as: abacavir, acyclovir, amprenavir, ascorbic acid vitamin C, azithromycin, baloxavir, boceprevir, cholecalciferol D, cidofovir, edoxudine, emtricitabine, and remdesivir. These ligands were analyzed using molecular docking, molecular quantum similarity, and chemical reactivity indices defined within a conceptual density functional theory framework. The analysis of molecular quantum similarity indices on inhibitors showed a high number of differences from a structural point of view. However, quite similar their electronic obtaining the highest values in the electronic similarity index. Global and local chemical reactivity stabilizing interactions using the crystal structure of SARS-CoV-2 RdRp. The molecular quantum similarity and chemical reactivity descriptors provide novel insights into these ligands that can be used in the design of new COVID-19 treatments.


Introduction
The severe acute respiratory syndrome coronavirus (SARS-CoV)-2 virus causes an infectious illness named coronavirus disease 2019 . SARS-CoV is a positive-sense single-stranded RNA virus from the Betacoronavirus genus. 1 Most RNA viruses require an RNA-dependent RNA polymerase (RdRp) for replication and transcription of the viral genome, thus making it essential for their survival, which is why it is a noteworthy topic to investigate. Its protein ranges from 240 to 450 kD and consists of a catalytic core and is a conserved protein within RNA viruses, and therefore it has been proposed as a potential option for the development of antiviral drugs. 2 Numerous drugs are being studied to treat COVID-19. In the USA, the first antiviral drug to treat COVID-19 in adults and children aged 12 years and older was remdesivir (Veklury). For hospitalized COVID-19 patients, remdesivir that is injected intravenously may be prescribed (Mayo Clinic).
Molecular docking, molecular quantum similarity (MQS), and global and local reactivity indices were used in this study to evaluate remdesivir and other related compounds such as abacavir, acyclovir, amprenavir, ascorbic acid vitamin C, azithromycin, baloxavir, boceprevir, cholecalciferol vitamin D, cidofovir, edoxudine, emtricitabine and hydroxychloroquine to obtain novel insights into how the process of stabilization of these ligands at the active site in the receptor structure takes place. Previous studies have demonstrated positive responses with different levels of effectiveness for these aforementioned compounds in the treatment of COVID-19 (Mayo Clinic).
Besalú et al., introduced MQS 30 years ago to study the similarity of molecules. 4 The density functional theory (DFT) framework is used to connect quantum mechanics and quantum chemistry utilizing molecular docking and chemical reactivity indices. 5 This research is leading to advances in the discovery of new treatment alternatives for COVID-19 and determining whether approved drugs, such as remdesivir, interact with other potential ligands.

Methodology for docking studies System preparation
The receptor structure, for the docking analysis, was extracted from the crystal structure of the SARS-CoV-2 RdRp (Protein Data Bank (PDB accession number, 6M71), which was adjusted using the Schrödinger (RRID:SCR_014879) suite 2017-1 Protein preparation Wizard (RRID:SCR_016749) module. i) The hydrogen bond network (H-bond) was optimized, and the protein structure refined, ii) at physiological pH, protonation states were determined using the PropKa utility, iii) The Impact Refinement (Impref) module was used to run a molecular minimization with heavy atoms constrained to a low root mean square deviation (RMSD) from the initial coordinates. [6][7][8] Moreover, the molecular structures of the compounds were sketched using Maestro Editor (Maestro, version 11.1, RRID: SCR_016748, Schrödinger, LLC). Then, 3D conformations were obtained with the LigPrep (RRID:SCR_016746) module, with ionization/tautomeric states predicted under physiological pH conditions with Epik (RRID:SCR_016745). Subsequently, energy minimization was used to comply with the Macro model using the OPLS2005 force field. The free alternative is AutoDock Vina (RRID:SCR_011958).

Molecular docking
Glide (RRID:SCR_000187) 9,10 with default parameters and Standard Precision (SP) model were used for docking investigations. The docking grid was created using default settings, with the co-crystallized ligand in the center. For the van der Waals radii of the nonpolar protein atoms, a scaling factor of 0.8 was applied to facilitate the binding of larger ligands. Extra precision (XP) was also utilized to expand alternate receptor conformations appropriate for binding ligands with unusual orientations via induced fit docking (IFD); this method allows the protein to undergo side-chain and/or backbone movements upon ligand docking. The results were optimized using the best predictions made by the 20 ns molecular dynamics.
Quantum similarity analysis Molecular quantum similarity measure A molecular quantum similarity measure (MQSM) amid two A and B systems, known as Z AB , compares two molecules using their respective density functions (DFs). Both DFs can be multiplied and integrated in terms of their electronic coordinates, which are then weighted using a predetermined positive operator Ω(r 1 ,r 2 ): [11][12][13] The nature of the operator used in Equation 1 provides the information being compared between the two systems and determines the similarity measure. For instance, if the chosen operator is the Dirac delta function (an efficient approach for functions with high peak values, such as the electronic density), i.e., Ω(r 1 ,r 2 )= δ(r 1 -r 2 ). One of the first similarity metrics employed is the overlapping MQSM; another widely used alternative is the Coulomb operator, i.e., Ω(r 1 ,r 2 )=| 1 -r 2 | -1 , resulting in a Coulombic MQSM. Even if the two molecules are equivalent, a similarity measure can be employed, which is known as a self-similarity measure (Z AA for the case of molecule A). 12 Given a group of N molecules, we can generate a similarity measure for each of them with regard to the other molecules in the group, including themselves. With this, a symmetric matrix can be constructed where the i-th column of the matrix is the collection of all similarity measures between the i-th molecule and each of the components of the group, including itself. Each vector (column of the matrix) is a discrete representation (in N dimensions) of the i-th structure. These sets of vectors are a set of chemical descriptors that do not simply express another set of molecular descriptors as is often done, but rather each descriptor has its own set of unique properties. 12,13 i. It is universal in that it can be derived from any collection of molecules and any individual molecule within that group.
ii. It is impartial because there are no other possibilities available in the construction process than those given the density functions and similarity measurements involved.

MQSM similarity index manipulation and visualization techniques
The similarity measure obtained for the group is unique once we have chosen a group of study objects and the operator in Equation 1; however, these measures can be combined or transformed to gain a new class of additional terms, which can be called Quantum Similarity Indices (QSI).
There are varieties of possible manipulations of the MQSM that result in a variety of QSI definitions. The most common ones are the following and are used in this paper [13][14][15][16] :

Carbó's similarity index
Carbó's similarity index between two molecules I and J is constructed from Equation 2. Because this index is also known as the cosine similarity index, it corresponds to the cosine of the angle formed by the density functions when considered as vectors. For any pair of compared molecules, this Carbo QSI has a value between 0 and 1, depending on the similarity between the two molecules (when I = J, the index approaches 1). [13][14][15][16][17] Quantum similarity using Euclidean distance Taking into account the similarity of Equation 3: It can be simplified to the so-called Euclidean distance index when k = x = 2 and can also be defined as follows: This Equation 4 forms the distance index of infinite order. 18

MQSM definition
Quantum similarity measurements are in accordance with psychological perception and the apparent principle of similarity: "the more similar two molecules, the more similar their properties are." This statement necessitates the construction of a molecule-to-molecule comparison, and the comparison of their densities is used to derive a quantitative measure of the quantum similarity of two molecules. Generally, the MQSM can be defined as the integral computational measure between two density functions, where the density functions are multiplied and integrated for the electronic coordinates in the relevant domain. [19][20][21] {r 1 , r 2 } are sets of electronic coordinates related to the corresponding wave function, and the operator Ω (r 1 ,r 2 ) is positively defined and supported on the electronic coordinates. [22][23][24][25][26] Types of measurements in molecular quantum similarity The types of measurements are mainly determined by the information required, the selection of supported operators and the form of MSQM. [18][19][20] MQSM overlap considering Equation 2: The distribution of Dirac's delta, Ω (r1, r2) = δ (r1, r2), is the most typical and intuitive choice for such a positively defined operator. These selections transform the broad definition of MQSM to compute the overlap MQSM that obtains measurements of the volume by both electronic density functions when they are superimposed. [17][18][19][20] The Dirac delta distribution, Ω (r1, r2) = δ (r1, r2), is the most typical and intuitive choice for a positive definite operator. This choice transforms the broad definition of MQSM to compute the overlap MQSM. [17][18][19][20] The MQSM calculates the degree of overlap between molecular densities using information about the electron "concentration" in the molecule [16][17][18][19][20][21] :

MQS coulomb considering Equation 2:
When the operator (Ω) is replaced with the Coulomb operator, Ω r 1 ,r 2 ð Þ¼ 1 jr1Àr2j , the coulomb MQS is generated, which defines the electrostatic Coulomb repulsion energy between two charge densities 20,21 : The Coulomb operator affects the overlap density functions. When considering molecular density functions as an electron distribution in space, this equation is simply an extension of the Coulomb operator for a distribution of continuous charge, thus can be used as electrostatic potential descriptors in some instances. This operator is correlated to electrostatic interactions and obtains the measurement of electrostatic repulsion between electronic distributions. [15][16][17][18][19] Euclidean distance index considering Equation 3: Another major transformation that can be expressed in terms of the classical distance is: Here Δx j ¼ x aj À x bj is the distance between a and b, and k = 2 is the definition of distance. Subsequently, the Euclidean distance between A and B (two quantum objects) is defined by: 17-21 Occasionally it is written as: where D AB has values in the range of [0,∞} but for situations where the compared items are identical, it converges to zero [17][18][19][20][21] : The norm of the differences of the density functions can be used in the definition Equation 9. The distance or dissimilarity index can be used to define the Euclidean distance index, which can also be represented as [15][16][17][18][19][20][21] : Alignment method: Topo-Geometrical Superposition Algorithm (TGSA) In this research, the TGSA 23 approach has been used to align the data. Gironés proposed the TGSA, in addition to programming and implementing it. This method assumes that the best way to align molecules is to superimpose them on a typical skeleton, with only the type atoms and the interatomic bonding interactions. The program accomplishes its purpose by examining the pairs of molecules and aligning the common substructure for a group. 23 Only topological and geometrical considerations are used in the procedure, including molecular topology and how distant bonds are compared. The superposition of two molecules is unique and is independent of the type of operator used to determine the similarity. 23,24 Chemical reactivity analysis Our research group has used several reactivity indices with Quantum Similarity. [26][27][28][29][30][31][32][33][34][35] The reactivity indices used in this work are chemical potential (μ), [36][37][38] hardness (ɳ), 39 and electrophilicity (ω), 38,39 which will be calculated. μ≈ The electrophilicity index (ω) measures the stabilization energy of the system when it is saturated by electrons from the external environment, and is mathematically defined as 38,39 : In this study, the local reactivity descriptors are the Fukui functions. Equations 15 and 16 represent the response of the chemical potential of a system to changes in the external potential. It is defined as the derivative of the electronic density with regard to the number of electrons at the constant external potential: Where f À k À Á is for nucleophilic attack and f À k À Á for the electrophilic attack. [40][41][42] In this approach, using the global and local reactivity descriptors makes it possible to compare the molecular reactivity at the sample set. All the structures were developed using M02X/6-31G(d, p) 43 methods in Gaussian 09 package. 44

Docking results
In this work, the stabilization process of a set of ligands related with activity against the novel SARS-CoV-2 was studied to understand the main ligand interactions in the active site, starting from the PDB crystal structure of SARS-CoV-2-dependent RNA polymerase code: 6M71. 45,46 The crystal structure of SARS-CoV-2 RdRp was selected taking into account the outcomes reported by Elfiky et al. 47 Figure 1 51 shows the docking interactions using remdesivir. In the conformation of Figure 1A and B, the docking interaction for the structure of remdesivir with a higher RMSD shows an -H bond with ARG555, ARG553 with two lengths of 2.28 Å and a length of 2.35 Å, LYS621 with a length of 2.22 Å, CYS622 with a length of 2.42 Å, and ASN691 with a length of 2.38 Å.
Another of the best conformations of remdesivir ( Figure 1C and D) involved a -H bond with residue ASN691 with a length of 2.46Å, two ARG553 and ARG555 with lengths of 2.35Å and 2.38Å, CYS622 with a length of 2.41Å, and LYS621 with a length of 2.39Å. The interactions with the residue ARG553 and ARG555 and LYS621 were very similar to Figure 1A and B. These interactions play an important role in the interaction of the active site, generating a good bonding surface (see Figure 1B and D). Figure 2 shows the docking interactions using ascorbic acid vitamin C structure of SARS-CoV-2 RdRp from SARS-CoV-2. The ligand and receptor surfaces show the active site zone defined by the stabilization interactions. A phosphate group was found to reside well inside a local binding pocket within the grip with residue CYS622.
In Figure 2A and B, the ascorbic acid vitamin C has -H bonds with the residues LYS621 at a distance of 2.20Å, ARG553 at a distance of 1.86Å, and LYS551 at a distance of 1.86Å, respectively. The ligand has the same interaction as remdesivir with residue LYS621. This conformation showed a good stabilization in the pocket, generating a good contour surface (see blue surface). Conformation 2 ( Figure 2C and D) shows three -H bonds with the same residue LYS621 at distances 2.18Å, 2.25Å and 2.27Å. This residue LYS621 is also shown in the conformations of remdesivir. Figure 3A and B shows the structure of cholecalciferol vitamin D, the largest RMSD corresponding to an -H bond with the residue SER759 at a length of 2.24Å. For conformation 2 in Figure 3C and D the -H bond is with the residue GLU166 with a length of 2.28Å. Finally, the compound has interactions with the residue TRP617 with a length 2.26 Å. Although this compound only interacted with the residue SER759, the ligand surface is well defined near the receptor, even though it is a single -H bond. The interaction is good enough to generate an active conformation with stabilizing capacity (Figure 3C and D), in agreement with previous reported results. 48 Conformation 1 of azithromycin ( Figure 4A and B), like those of remdesivir and ascorbic acid vitamin C, has -H bond with LYS621 with a length of 2.28Å. In addition, it has two interactions with the residue ASP760 at lengths of 2.26Å and 2.34Å, with ASP761 at a length of 2.32Å. On the other hand, for conformation 2 of azithromycin of Figure 4C and D, we can see two -H bonds with the residues ASP760 with lengths of 2.41Å and 2.38Å, ASP761 with a length of 2.38Å and TRP617 with a length of 2.28Å. These ASP760, ASP761 and TRP617 interactions agree with the recent work published by Abdrabbo et al. 49 In this work, once the stable binding mode was located, free energy perturbation (FEP) calculations were carried out to estimate the binding affinity of RemTP and ATP to COVID-19, and to identify the key residues in the binding process.

Conformation 1 of
The docking interactions in conformation 1 of Figure 5A and B for hydroxychloroquine with higher RMSD shows a -H bond with the residue ARG553 with a length of 2.38Å, with ASP760 at lengths of 2.36 Å, 2.28 Å and 2.32 Å. For the best  Figure 5C and D involved ASP462, ASP623 and ASP760 with lengths of 2.36 Å, 2.31 Å and 2.28 Å, respectively.
In this study, a molecular quantum similarity study was developed to systematically investigate and analyze the structural and electronic considerations involved in active site stabilization. An electronic similarity study has been designed to find the relationship between ligands and the active site of RdRp proteins.

Molecular quantum similarity indices
Because conformations are crucial in docking studies, we have analyzed the impact that molecular alignment has on the generation of interactions (-H bond) at the active site. These molecular alignments, which generated the best poses, were related to potential COVID-19 inhibitors. For this analysis we have used four similarity descriptors (structural, electronic and their Euclidean distances) and molecular quantum similarity indices were calculated using the overlap and coulomb descriptors (Equations 6 and 7).
The compounds analyzed in the docking results were compared with a series of ligands with SARS-CoV-2 activity such as abacavir, acyclovir, amprenavir, ascorbic acid vitamin C, azithromycin, baloxavir, boceprevir, cholecalciferol vitamin D, cidofovir, edoxudine, emtricitabine, hydroxychloroquine and remdesivir. These compounds have been associated with COVID-19 treatment.
The molecular quantum similarity indices are defined on the interval (0,1), where 0 means dis-similarity and 1 means selfsimilarity. Table 1 shows the molecular quantum similarity indices using the overlap descriptor. Based on the distribution of Dirac's delta, Ω (r1, r2) = δ (r1, r2), this operator can be related to the structural characteristics of the compound analyzed. The highest values were for abacavir/acyclovir (0.6286) with an Euclidean distance of 3.5612 (Table 2)  The molecular quantum similarity indices using the Coulomb descriptor were calculated to obtain insights into the electronic similarity (see Table 3). The Coulomb descriptor is related to the electronic similarity in the molecular set. Table 3 shows that the highest values for emtricitabine/edoxudine (0.9390) with an Euclidean distance of 14.4795 (see Table 4), emtricitabine/abacavir (0.9311) with an Euclidean distance of (16.2459) and edoxudine/acyclovir (0.9165) with a Euclidean distance of 22.2556. Unlike the overlap descriptor, the values are near to 1. This means that despite having significant structural differences, these structures have high electronic similarities.

Global and local reactivity indices
This work has also explored the global and local chemical reactivity indices within the DFT framework. Table 5 shows that the least reactive molecule was hydroxychloroquine, whose electronic chemical potential μ= -1.7070 eV, chemical hardness Ƞ=1.606 eV, softness S=0.6226 eV and electrophilicity ω=0.4536 eV. The most reactive compound was baloxavir-marboxil hydroxychloroquine, whose electronic chemical potential μ= -3.7335 eV, chemical hardness Ƞ=2.1441 eV, softness S=0.4664 eV and electrophilicity ω=1.6253 eV. These electrophilicity values play an essential role in ligands stabilized by non-covalent interactions and determine the stability of the active site.

1.0000
O_Hab, overlap index; remd, remdesivir; abaca, abacavir; acyc, acyclovir; ampre, amprenavir; asco, ascorbic acid vitamin C; azythr, azithromycin; balox, baloxavir; boce, boceprevir; chole, cholecalciferol vitamin D; cido, cidofovir; edox, edoxudine; emitri, emtricitabine; hidrox, hydroxychloroquine.    Figure 6 shows the Fukui functions that have been used to describe the local chemical reactivity of remdesivir and cholecalciferol vitamin D. Figure 6 and Figure 7 show the Fukui functions f À r ð Þ≈ HOMO j j 2 and f þ r ð Þ≈ HOMO j j 2 for the compounds selected. The Fukui function f À r ð Þ≈ HOMO j j 2 is associated with the susceptibility of a site to being attacked by electrophilic species, while f þ r ð Þ≈ HOMO j j 2 represents the susceptibility of a site to being attacked by nucleophilic species. In these figures, we can see the functions of Fukui (-) and (+) on the identical zones in the molecular structure. These can be related to the stabilization process on the active site, showing the zones associated with global reactivity indices like electrophilicity and softness. These reactivity parameters can help understand the destabilization process and the bond formation (-H bonds) in non-covalent interactions. All these outcomes can be useful for the rapid assessment of the currently available antiviral drugs used for treating COVID-19 patients, which is crucial at this time of crisis, as well as for discovering newer drugs.

Conclusions
In this study, a number of ligands related to the treatment of COVID-19 and used as inhibitors of SARS-CoV-2 virus, such as abacavir, acyclovir, amprenavir, ascorbic acid vitamin C, azithromycin, baloxavir, boceprevir, cholecalciferol vitamin D, cidofovir, edoxudine, emtricitabine, hydroxychloroquine and remdesivir, have been analyzed by molecular docking, molecular quantum similarity and chemical reactivity indices to study their active site stabilization interactions from a structural and electronic point of view.
In this study, to develop the docking experiments, some compounds were selected as references, these compounds are frequently related with the treatment of COVID-19. Some important anti-polymerase drugs were tested that are currently in clinical trials or on the market to stop viral infection on an emergency basis. The docking interaction for remdesivir, cholecalciferol vitamin D, azithromycin and ascorbic acid have shown good interaction (-H bonds) in the active site. The main idea is to extrapolate these outcomes and find novel insights into inhibitors for COVID-19.
The analysis of molecular quantum similarity indices on inhibitors showed high differences from a structural point of view. However, they are quite similar in their electronic density, obtaining the highest values in the electronic similarity index. Global and local chemical reactivity indices were calculated. These indices allow for the identification of the zones of chemical reactivity involved in the non-covalent stabilization of these inhibitors on the active site. Moreover, new outcomes about compounds such as abacavir, which is used to treat HIV infection, were shown and some studies have shown their relationship with COVID-19. In addition, it sheds light on the use of novel ligands for the treatment of this challenging disease that has claimed millions of lives worldwide.

Data availability
Underlying data Protein Data Bank: SARS-Cov-2 RNA-dependent RNA polymerase in complex with cofactors. Accession number 6M71; https://identifiers.org/pdb:6M71 Harvard Dataverse: New insights in series of ligands used as SARS-CoV-2 virus inhibitors within molecular docking, molecular quantum similarity, and chemical reactivity indices frameworks. https://doi.org/10.7910/DVN/MSIGUS. 51 Data are available under the terms of the Creative Commons Zero D. "No rights reserved" data waiver (CC0 1.0 Public domain dedication).