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

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–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₁,r₂):^11–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₁,r₂)= δ(r₁- r₂). One of the first similarity metrics employed is the overlapping MQSM; another widely used alternative is the Coulomb operator, i.e., Ω(r₁,r₂)=|‌₁- r₂|‌‌^-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).¹²

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

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–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–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.¹⁸

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–21

{r₁, r₂} are sets of electronic coordinates related to the corresponding wave function, and the operator Ω (r₁,r₂) is positively defined and supported on the electronic coordinates.^22–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–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–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–20 The MQSM calculates the degree of overlap between molecular densities using information about the electron “concentration” in the molecule^16–21:

MQS coulomb considering Equation 2:

When the operator (Ω) is replaced with the Coulomb operator, $\mathrm{\Omega }\left({\mathrm{r}}_1,{\mathrm{r}}_2\right)=\frac{1}{|r_1-r_2|}$, 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–19

Euclidean distance index considering Equation 3:

Another major transformation that can be expressed in terms of the classical distance is:

Here $\mathrm{\Delta }{x}_j=x_{\mathit{aj}}-x_{\mathit{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: $D_{\mathit{AB}}=\sqrt{Z_{\mathit{AA}}+Z_{\mathit{BB}}+Z_{\mathit{AB}}}$, where D_AB has values in the range of [0,∞} but for situations where the compared items are identical, it converges to zero^17–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–21:

Alignment method: Topo-Geometrical Superposition Algorithm (TGSA)

In this research, the TGSA²³ 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.²³ 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

Molecular alignment is a critical factor in molecular quantum similarity studies, which aim to compare and analyze the similarity of molecules based on their quantum mechanical properties. Molecular alignment refers to the arrangement or orientation of molecules in a specific way that facilitates the comparison of their electronic structure and other quantum properties. The importance of molecular alignment in molecular quantum similarity can be understood through the following points: Consistent Geometries: Molecular alignment ensures that the molecules being compared are in consistent geometries or orientations. This is crucial because the electronic structure and quantum properties of molecules can vary significantly based on their spatial arrangement. By aligning molecules in a specific way, researchers can focus on comparing their intrinsic quantum properties rather than being confounded by geometric differences. Accurate Comparison: Quantum similarity aims to quantify the resemblance of molecular properties, such as electronic distributions and energy levels. Proper alignment enables a more accurate comparison of these properties across different molecules. It allows researchers to assess whether molecules with similar alignments also exhibit similar quantum behavior, providing insights into their functional and structural relationships. Quantitative Analysis: Molecular quantum similarity often involves complex computational methods to calculate molecular properties like electronic densities, molecular orbitals, and energy levels. Accurate alignment ensures that these calculations are based on consistent molecular arrangements, leading to reliable quantitative analyses of similarity. Drug Discovery and Materials Design: Molecular quantum similarity is employed in drug discovery and materials design to identify compounds with similar quantum properties to known active compounds or materials. Proper molecular alignment ensures that the comparison is meaningful, leading to the identification of candidates with potential similar functional properties. Structure-Property Relationships: Understanding how molecular alignment affects quantum similarity can provide insights into the relationship between molecular structure and properties. Researchers can gain a deeper understanding of how specific arrangements influence quantum properties and functional behavior. Enhancing Predictive Models: Machine learning and computational methods often rely on molecular quantum similarity to predict various properties of molecules. Accurate molecular alignment contributes to the quality of training data and helps in building more reliable predictive models. Quantum Chemistry Interpretations: Proper molecular alignment facilitates the interpretation of quantum chemical calculations. It helps researchers visualize and understand the impact of different molecular orientations on electronic structure, energy levels, and other quantum properties.

In summary, molecular alignment is crucial in molecular quantum similarity studies because it allows for consistent and accurate comparisons of quantum properties across different molecules. This alignment ensures that the observed similarities are not confounded by variations in molecular orientation, leading to meaningful insights into the functional and structural relationships between molecules.

Chemical reactivity analysis

Our research group has used several reactivity indices with Quantum Similarity.^26–35 The reactivity indices used in this work are chemical potential (μ),^36–38 hardness (ɳ),³⁹ 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 $\left(f_k^{+}\right)$ is for nucleophilic attack and $\left(f_k^{-}\right)$ for the electrophilic attack.^40–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)⁴³ methods in Gaussian 09 package.⁴⁴

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.⁴⁷ Figure 1⁵¹ 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 Å.

Figure 1. Docking outcomes using remdesivir.

A. Docking interactions using remdesivir with the crystal structure of SARS-CoV-2 RdRp (PDB accession number 6M71). B. Docking interactions using Receptor (gray) and ligand (blue) surfaces. C. Docking interactions with the receptor site. D. Surface of the binding pocket with the receptor site. SARS-CoV-2, severe acute respiratory syndrome coronavirus 2; RdRp, RNA-dependent RNA polymerase; PDB, Protein Data Bank.

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.

Figure 2. Docking outcomes using ascorbic acid vitamin C.

A. Docking interactions using ascorbic acid vitamin C with the crystal structure of SARS-CoV-2 RdRp (PDB accession number 6M71). B. Receptor (gray) and ligand (blue) surfaces. C. Docking interactions with the receptor site. D. Surface of the binding pocket with the receptor site. SARS-CoV-2, severe acute respiratory syndrome coronavirus 2; RdRp, RNA-dependent RNA polymerase; PDB, Protein Data Bank.

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.

Conformation 1 of 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.⁴⁸

Figure 3. Docking outcomes using cholecalciferol vitamin D.

A. Docking interactions of cholecalciferol vitamin D with the crystal structure of SARS-CoV-2 RdRp (PDB accession number 6M71). B. Receptor (gray) and ligand (blue) surfaces. C. Docking interactions with the receptor site. D. Surface of the binding pocket with the receptor site. SARS-CoV-2, severe acute respiratory syndrome coronavirus 2; RdRp, RNA-dependent RNA polymerase; PDB, Protein Data Bank.

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.⁴⁹ 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.

Figure 4. Docking outcomes using azithromycin.

A. Docking interactions using azithromycin with the crystal structure of SARS-CoV-2 (PDB accession number 6M71). B. Receptor (gray) and ligand (blue) surfaces with a receptor (gray) and ligand (blue) surfaces. C. Docking interactions with the receptor site. D. Surface of the binding pocket with the receptor site. SARS-CoV-2, severe acute respiratory syndrome coronavirus 2; RdRp, RNA-dependent RNA polymerase; PDB, Protein Data Bank.

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 conformation 2, Figure 5C and D involved ASP462, ASP623 and ASP760 with lengths of 2.36 Å, 2.31 Å and 2.28 Å, respectively.

Figure 5. Docking outcomes using hydroxychloroquine.

A. Docking interactions using hydroxychloroquine with the crystal structure of SARS-CoV-2 RdRp (PDB accession number 6M71). B. Receptor (gray) and ligand (blue) surfaces with a receptor (gray) and ligand (blue) surfaces. C. Docking interactions with the receptor site. D. Surface of the binding pocket with the receptor site. SARS-CoV-2, severe acute respiratory syndrome coronavirus 2; RdRp, RNA-dependent RNA polymerase; PDB, Protein Data Bank.

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 self-similarity. 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), emtricitabine/abacavir (0.6011) with an Euclidean distance of 3.7645 and emtricitabine/edoxudine (0.6079) with a Euclidean distance of 3.8020. The lowest values were for boceprevir/azithromycin (0.1374), and remdesivir/cholecalciferol (0.1476) with a Euclidean distance of 7.1944.

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.

The compounds with the highest chemical potential (negative electronegativity) were hydroxychloroquine with μ=-1.7070 eV, abacavir with μ=-2.6330 eV and cholecalciferol vitamin D with μ=-2.9825 eV. These compounds are commonly used in COVID-19 treatment.^48–50 On the other hand, abacavir is used to treat HIV infection and some studies have shown its relationship with COVID-19.⁵⁰ However, some countries are skeptical about its use for COVID-19 treatment. Other important compounds often used for the treatment of SARS-CoV-2 virus are azithromycin and remdesivir, with a chemical reactivity of μ=-3.5621 eV and -3.6097 eV, chemical hardness η=2.5532 eV and 2.3253 eV, softness S=0.3916 eV and 0.4300 eV and electrophilicity ω=1.2136 eV and 1.4009 eV, respectively. Both are very reactive and have high electrophilicity values. Figure 6 shows the Fukui functions that have been used to describe the local chemical reactivity of remdesivir and cholecalciferol vitamin D.

Figure 6. Local chemistry reactivity indices.

A. Fukui function $f^{-}\left(r\right)\approx {\left|\mathit{\text{HOMO}}\right|}^2$ and B. Fukui Function $f^{+}\left(r\right)\approx {\left|\mathit{\text{HOMO}}\right|}^2$ related with remdesivir. C. Fukui Function $f^{-}\left(r\right)\approx {\left|\mathit{\text{HOMO}}\right|}^2$ and D. $f^{+}\left(r\right)\approx {\left|\mathit{\text{HOMO}}\right|}^2$ related with cholecalciferol vitamin D.

Figure 6 and Figure 7 show the Fukui functions $f^{-}\left(r\right)\approx {\left|\mathit{\text{HOMO}}\right|}^2$ and $f^{+}\left(r\right)\approx {\left|\mathit{\text{HOMO}}\right|}^2$ for the compounds selected. The Fukui function $f^{-}\left(r\right)\approx {\left|\mathit{\text{HOMO}}\right|}^2$ is associated with the susceptibility of a site to being attacked by electrophilic species, while $f^{+}\left(r\right)\approx {\left|\mathit{\text{HOMO}}\right|}^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.

Figure 7. Local chemistry reactivity indices.

A. Fukui function $f^{-}\left(r\right)\approx {\left|\mathit{\text{HOMO}}\right|}^2$ and B. $f^{+}\left(r\right)\approx {\left|\mathit{\text{HOMO}}\right|}^2$ related with ascorbic acid vitamin C. C. Fukui function $f^{-}\left(r\right)\approx {\left|\mathit{\text{HOMO}}\right|}^2$ and D. $f^{+}\left(r\right)\approx {\left|\mathit{\text{HOMO}}\right|}^2$ related with azithromycin. E. Fukui function $f^{-}\left(r\right)\approx {\left|\mathit{\text{HOMO}}\right|}^2$ and F. $f^{+}\left(r\right)\approx {\left|\mathit{\text{HOMO}}\right|}^2$ associated with hydroxychloroquine.

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.⁵¹

Data are available under the terms of the Creative Commons Zero D. “No rights reserved” data waiver (CC0 1.0 Public domain dedication).