Study anti-viral drugs for their efficiency against multiple SARS CoV-2 drug targets within molecular docking, molecular quantum similarity, and chemical reactivity indices frameworks

2.1 System preparation

In the docking experiment, the receptor structure was derived following specific protocols based on the crystal structure of SARS-CoV-2 RNA-dependent RNA polymerase with PDB code 6m71. Adjustments to the structure were made using the protein preparation wizard module from the Schrödinger suite 2017-1. These adjustments included:

Conversely, the molecular structures of the compounds were constructed using the Maestro Editor (Maestro, version 11.1, Schrödinger, LLC). Subsequently, 3D conformations were generated using the LigPrep module, and ionization/tautomeric states were predicted under physiological pH conditions using Epik. Finally, energy minimization was performed with the Macro model using the OPLS2005 force field.

2.2 Molecular docking

Glide (https://newsite.schrodinger.com/platform/products/glide/)^9,10 with default parameters and Standard Precision (SP) model was 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, backbone, or both movements upon ligand docking. All results were redocking and RMSD were performed. The binding pocket of the RdRp—GLY616, TRP617, ASP618, TYR619, LEU758, SER759, ASP760, ASP761, ALA762, LYS621, TYS799, TRP800, GLU811, PHE812, CIS813, and/or SER814—was found using Glide.^9,10

The docking process involves four precise steps, relying on Glide’s scoring function and Prime’s advanced conformational refinement to ensure accuracy:

To further assess the interactions of the ligands in the active site, the extent of residue movement induced by the IFD computation is considered. For both the most and least active ligands, all poses are compared within the molecular set. Molecular dynamics calculations over 30ns are employed to predict the best poses and analyze their stabilization in the active site.

3.1 Molecular Quantum Similarity Measure

A Molecular Quantum Similarity Measure (MQSM) amid two A and B systems, known as Z_AB, compares two molecules that may be created 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 operator used in Equation 1 plays a crucial role in determining the information being compared and serves as the measure of similarity between the two systems. For example, when the operator chosen is the Dirac delta function, it proves to be an efficient approach for functions with high peak values, like the electronic density. Moreover, it provides a similar mathematical abstraction as a charge or point mass, 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₁ - r₂|‌‌⁻¹, resulting in a Coulombic MQSM. Even if the two molecules are equivalent, a similarity measure can be employed for any two molecular systems; this measurement is known as a self-similarity measure (Z _AA for the case of molecule A).¹²

For a given group of N molecules, we can derive a measure of similarity for each molecule concerning the others in the group, including itself. These similarity measurements can then be used to construct a symmetric matrix. The i-th column of this matrix represents a compilation of all similarity measurements between the i-th molecule and every constituent in the group, including itself. Consequently, each vector (matrix column) serves as a discrete N-dimensional representation of the i-th structure. These vector sets can be demonstrated as a set of chemical descriptors. However, this set of similarity matrix columns goes beyond merely representing another set of molecular descriptors, as commonly done for theoretical molecule description; each descriptor possesses unique properties.^12–20

3.2 Carbó’s similarity index

Carbó’s similarity index between two molecules I and J are constructed from Equation 2. Because this index is also known as the cosine similarity index, it corresponds to the cosine of the angle included by the density functions involved 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–28

3.3 The quantum similarity using the Euclidean distance

Taking into account the similarity of equation 3:

It is simplified to the so-called Euclidean distance index when k = x = 2. Index 3 of the form can also be defined as follows:

This Equation 4 forms the distance index of infinite order.^18–30

3.4 MQSM overlap considering the Equation 2

The distribution of Dirac’s delta, Ω (r₁, r₂) = δ (r₁, r₂), 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 function is derived instinctively from the physical definition, and it is computationally compliant. The MQSM calculates the degree of overlap between molecular comparisons using information about the electron concentration in the molecule.^16–21

3.5 MQS Coulomb considering the Equation 5

When the operator (Ω) is replaced with the Coulomb operator, Ω (r₁, r₂)= $\frac{1}{│r_1-r_2│}$ , the coulomb MQS is generated, which defines the electrostatically repellent coulomb 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 the 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.^30–37

3.6 Euclidean distance index considering the Equation 3

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

Here $\Delta x_j=x_{\mathit{aj}}-x_{\mathit{bj}}$ is the distance between a and b, and k = 2 is the definition of distance, respectively. Subsequently, the Euclidean distance between A and B as two quantum objects are defined by^17–21:

Occasionally it is written as: $D_{\mathit{AB}}=\sqrt{Z_{\mathit{AA}}+Z_{\mathit{BB}}+{\mathit{ZZ}}_{\mathit{AB}}}$ , where D _AB has values in the range of [0,∞) but for earlier circumstances where the compared items are identical, it converges to zero between them.^17–21:

The norm of the differences in the density functions of the compared objects can be used to interpret this index geometrically. The distance or dissimilarity index can be used to define the Euclidean distance index, which can also be represented as^21–25:

3.7 Alignment method: Topo-Geometrical Superposition Algorithm (TGSA)

In this investigation, the TGSA (Typical Geometry Superposition Algorithm) approach was utilized for data alignment. Devised by Gironés, TGSA operates under the assumption that the optimal way to align molecules involves superimposing them onto a shared skeleton, considering solely the atomic types and interatomic bonding interactions based on the atomic number coordination.^23–28

The program initiates by scrutinizing pairs of atoms in the molecules, aligning their common substructure for a group of molecules using topological and geometrical considerations. Notably, the superposition achieved is distinctive and unaffected by the choice of similarity measure.^28–34

Initially, the program organizes molecular coordination into bases based on the reduction of atomic numbers, defining a path for the number of hydrogens in the molecule (excluding hydrogen atoms for computational efficiency). Subsequently, atomic pairs are delineated, specifying the involved atoms and their respective distances.^34–37 Translocations are identified through changes in the conformations’ spine caused by substitutions in the molecules.²³ Bones not aligning with the skeletons are eliminated during this process.

The program then assembles atomic triads by incorporating three atoms from the compared pairs, forming a triangle in the plane representing the chemical box’s efficacy.^23,24 The triangles generated for two molecules are compared using their respective interatomic and translational distances. Triads meeting the classification criteria are retained, superimposed, and dictate the molecular alignment result.²³

It’s worth noting that since TGSA characterizes molecules as rigid structures without flexibility (no vibration or rotation in bond distances and angles), it may not yield optimal results for diverse molecular structures due to the restricted alignment with the common recognition skeleton. Nevertheless, this method consistently aligns with chemical intuition and is favored for its accessibility and lower computational requirements.^26–37

MQS can help in optimizing the drug-likeness of a molecule by evaluating quantum-level descriptors that influence the bioavailability, stability, and efficacy of drugs. For example, electron distribution affects how a molecule will interact with biological macromolecules, and using MQS can help modify structures to improve these interactions.

5.1 Molecular docking outcomes for the ligands

Drug discovery for COVID-19 presents several challenges, and researchers worldwide have been working diligently to address these issues. Some of the key challenges include: Virus Mutability: SARS-CoV-2, the virus responsible for COVID-19, can mutate, leading to the emergence of new variants. This mutability poses a challenge in developing drugs that can effectively target different strains of the virus. Rapid Evolution of the Pandemic: The rapid spread of the virus and the urgent need for effective treatments make it challenging to follow traditional drug development timelines. Accelerated timelines can compromise thorough testing and validation processes. Lack of Pre-existing Therapies: Unlike some other infectious diseases, there were no pre-existing drugs specifically designed to target SARS-CoV-2. Developing new drugs from scratch is a time-consuming process.

Complexity of the Virus Life Cycle: Understanding the intricate details of the virus’s life cycle and the host-pathogen interactions is essential for developing targeted therapies. This complexity requires a deep understanding of virology and immunology. Drug Safety: Ensuring the safety of potential treatments is crucial. Some drugs may show promise in early stages but could have adverse effects or interactions with other medications, requiring extensive testing for safety. Drug Delivery Challenges: Designing effective drug delivery systems to ensure that the drug reaches the target tissues in sufficient concentrations is a significant challenge. This is especially important for antiviral drugs targeting the respiratory system. Antibody Resistance: The virus may develop resistance to certain antiviral drugs or antibodies over time. This highlights the need for ongoing research to identify multiple targets for drug development and combination therapies. Global Collaboration: International collaboration is crucial for sharing data, resources, and expertise. However, coordinating efforts across borders and overcoming logistical and political challenges can be complex. Vaccine Success and Impact: The success and widespread distribution of COVID-19 vaccines have been crucial in controlling the pandemic. However, ongoing research is needed to address vaccine effectiveness against new variants and to develop treatments for breakthrough cases. Resource Limitations: Drug discovery requires significant financial and human resources. The COVID-19 pandemic has strained healthcare systems globally, and prioritizing and allocating resources for research and development can be challenging. Despite these challenges, the scientific community has made remarkable progress in a short time, developing vaccines and exploring various therapeutic approaches. Continuous research and collaboration will be essential for addressing the evolving nature of the pandemic and improving our ability to manage and treat COVID-19. For these reasons, in this study are obtained new insights for a serie of ligands, based on the crystal structure of SARS-CoV-2 RNA-dependent RNA polymerase with PDB code 6m71. Please, see Figures 1 and 2.

Figure 1. Docking results for Oseltamivir.

Figure 2. Docking results for Prochloraz.

In Figure 1 are shows the docking outcomes for Oseltamivir. The main interaction is with the residues GLU166 (-H, 1.323 Å), GLN189 (-H, 1.185 Å) and GLY143 (-H, 1.054 Å). However, in Figure 2 for Prochloraz has a stacking (also called π–π stacking, (1.213 Å)) with the residues HIE41.

The Figure 3 shows the docking results for Valacyclovir, this ligand shows interactions with the residues CYS145 (-H, 1.452 Å), SER144 (-H, 1.114 Å), GLY143 (-H, 1.156 Å) and HIS163 (-H, 1.254 Å).

Figure 3. Docking results for Valacyclovir.

In Figure 4 show the interaction for Baricitinib, this Figure shows interactions with the residues GLU166 (-H, 1.568 Å), GLY143 (-H, 1.365 Å). Unlike, in the Figure 5 we can see interactions for Molnupiravir with the residue CYS145 (-H, 1.248 Å and 1.238 Å).

Figure 4. Docking results for Baricitinib.

Figure 5. Docking results for Molnupiravir.

In Figure 6, shows the interactions for Valacyclovir with the residue HIS163 (-H, 1.156 Å) and CYS145 (-H, 1.256 Å). The Figure 7, shows the interactions for Penciclovir with the residues HIS164 (-H, 1.005 Å) and π-π interactions with the residues HIE41 (-H, 1.269 Å).

Figure 6. Docking results for Valacyclovir.

Figure 7. Docking results for Penciclovir pose 1.

In the Figures 8 and 9, shows the interactions with the residues HIE41 (-H, 1.567 Å) for Penciclovir and with the residues GLY143 (-H, 1.456 Å), CYS145 (-H, 1.436 Å) and HIS163 (-H, 1.485 Å).

Figure 8. Docking results for Penciclovir pose 2.

Figure 9. Docking results for Famciclovir.

Finally, in the Figures 10 and 11 we can see the interactions with the residues GLU166 (-H, 1.054 Å) and GLY143 (-H, 1.158 Å) for Lamivudine and with the residues GLU166 (-H, 1.254 Å), HIE41 (-H, 1.266 Å) and GLY143 (-H, 1.354 Å) for Nitazoxanide.

Figure 10. Docking results for Lamivudine.

Figure 11. Docking results for Nitazoxanide.

5.2 Molecular quantum similarity analysis

MQSM helps in the identification of molecules with similar quantum properties, which can be useful for finding potential drug candidates. Similar quantum properties suggest similar chemical behaviours. MQSM involves the calculation of quantum descriptors, which are numerical representations of molecular properties derived from quantum mechanical calculations. These descriptors may include electronic density, electron localization function, and others. Various similarity measures, such as Carbo’s similarity index or Euclidean distance, are employed to quantify the degree of similarity between quantum descriptors of different molecules MQSM can aid in lead optimization by identifying compounds with quantum properties like known drug candidates. This assists in designing molecules with improved pharmacological profiles.

The MQSM can be used to understand the interactions between ligands and their target receptors at a quantum level. This is crucial for rational drug design. It is often used in conjunction with molecular docking studies. While docking predicts the binding affinity and geometry of ligands with target proteins, MQSM provides insights into the quantum properties that influence these interactions.

MQSM offers a detailed and atomistic understanding of molecular properties, allowing for a more nuanced analysis of chemical similarity. Unlike traditional structural similarity measures, MQSM considers the quantum properties of molecules, providing a more comprehensive comparison.

In the context of antiviral drug discovery, MQSM can be employed to identify molecules with similar quantum properties to known antiviral drugs, aiding in the search for new therapeutics. Therefore, Molecular Quantum Similarity analysis is a valuable tool in drug discovery that leverages quantum mechanical principles to assess the similarity of molecular properties. This approach contributes to the rational design of novel drugs and the optimization of lead compounds for improved pharmacological profiles.

In Tables 1 and 2. the higher MQSM is between the compounds Lamivudine and Molnupiravir using the Overlap operator 0,5742 with a euclidean distance of 4,2364, see Table 2. The lowest MQSM is between the compounds Oseltamivir and Prochloraz (0,2233) with a eclidean distance of 5,5841 (see Table 2). These measurements are below 0.5, according to the range of carbon indices (0.1]. So we can say that structurally they are quite different.

In Tables 3 and 4 the higher MQSM using the Coulomb operator is between the compounds Lamivudine and Molnupiravir (0,9178) with a euclidean distance of 22,5461, see Table 4. However, the lowest MQSM is between the compounds Baricitinib and Farmciclovir (0,6001) with a euclidean distance of 44,3298. The comparison between the compounds Lamivudine and Molnupiravir has the higher values in both case using the Overlap and Coulomb operators. Unlike to the MQSM using the overlap operator, the measurements of the Table 3 are above 0.5, according to the range of carbon indices (0.1]. Therefore, we can say that although structurally they are not so similar, electronically they are. Because the electronic similarity index is higher with respect to the overlap, we have analyzed in depth the electronic properties of the analyzed ligands.

5.3 Global reactivity descriptors analysis and Fukui function comparison

The investigation explored into global and local chemical reactivity descriptors through DFT computations. This part contrasts the reactivity of the ligands within the study, encompassing both overarching parameters and locally descriptive functions of reactivity. Electrophilicity values hold potential significance in stabilizing the active site of ligands engaged in non-covalent interactions. Illustrated in Figure 12 are the computed global parameters—chemical potential, chemical hardness, global softness, and global electrophilicity—offering a comparative analysis of the ligand sample’s chemical reactivity. As depicted in Figure 12, Nitazoxanide emerges as the most reactive molecule, displaying the lowest values of electronic chemical potential (μ) and chemical hardness (η), alongside the highest global softness and global electrophilicity. Conversely, the remaining compounds exhibit descriptor values that display less pronounced differences. Consideration of solely global descriptors might suggest comparable reactivity. The values obtained in each case are detailed in Table S1 in the accompanying supplementary information.

Figure 12. Global parameters, including electronic chemical potential (A), chemical hardness (B), global softness (C), and global electrophilicity (D) for the compounds in the study.

Since the analysis of the global parameters is limited, we will complete it with the comparison of some local descriptor functions. The electrophile and nucleophile Fukui functions (as a measure of reactivity) were then compared using the Frontier Molecular Orbital (FMO) approach. The electrophilic-nucleophilic character of the following functions also shows those molecular areas that are most likely to form charge-donating interactions (basically by charge delocalisation). These types of interactions are important and difficult to determine using docking analysis. Figure 13 shows Fukui function $f^{-}\left(\overrightarrow{r}\right)$ calculated under the FMO approximation ( ${\left|\mathit{\text{HOMO}}\left(\overrightarrow{r}\right)\right|}^2$ ) for the compounds Valacyclovir and Penciclovir (A and B respectively), it can be noted that $f^{-}\left(\overrightarrow{r}\right)$ is similar in both compounds, and in both cases, the function assigns the most nucleophilic character to the condensed rings. In Figures S1-S8, see Underlying Data, the $f^{-}\left(\overrightarrow{r}\right)$ functions for all ligands in the study can be seen. Comparing the ${\left|\mathit{\text{HOMO}}\left(\overrightarrow{r}\right)\right|}^2$ function from Figure 13-B with Figure 8, it can be observed that the interactions of Penciclovir with GLY 143, HIE 41, and LEU 141 obtained through docking are compatible with the $f^{-}\left(\overrightarrow{r}\right)$ function.

Figure 13. Fukui function $f^{-}\left(\overrightarrow{r}\right)$ calculated under the FMO approximation $\left({\left|\mathit{\text{HOMO}}\left(\overrightarrow{r}\right)\right|}^2\right)$ for the ligands A) Valacyclovir and B) Penciclovir.

Isovalue was 0.008 in both cases. The figure was created using GaussView 5.0.

Figure 14 the functions $f^{+}\left(\overrightarrow{r}\right)$ calculated under the FMO approximation ( ${\left|\mathit{\text{LUMO}}\left(\overrightarrow{r}\right)\right|}^2$ ) for compounds A) Prochloraz, B) Molnupiravir and C) Lamivudine. It can be noted that $f^{+}\left(\overrightarrow{r}\right)$ is similar in the three compounds; in all three cases, the function assigns the most nucleophilic character to various carbons in the six-carbon ring. The ring may have some ability to rotate and orient itself towards where it can form the strongest interactions. Comparing Figure 14 with Figure 2, it can be seen that $f^{+}\left(\overrightarrow{r}\right)$ is compatible with the interaction of Prochloraz with HIE 41. A similar conclusion is drawn when comparing Figure 14 with Figure 5, in this case, the interaction is formed between Molnupiravir and PHE 140. Lastly, comparing Figure 14 with Figure 10 highlights the interaction between Lamivudine and GLU 166.

Figure 14. Fukui function $f^{+}\left(\overrightarrow{r}\right)$ calculated under the FMO approximation $\left({\left|\mathit{\text{LUMO}}\left(\overrightarrow{r}\right)\right|}^2\right)$ for the ligands A) Prochloraz, B) Molnupiravir and C) Lamivudine.

The isovalue was 0.008 in all cases. The figure was created using GaussView 5.0.

Figure 15 shows the functions $f^{+}\left(\overrightarrow{r}\right)$ calculated under the FMO approximation ( ${\left|\mathit{\text{LUMO}}\left(\overrightarrow{r}\right)\right|}^2$ ) for compounds A) Valacyclovir and B) Penciclovir. It can be noted that $f^{+}\left(\overrightarrow{r}\right)$ is similar in the two compounds, and in both cases, the function assigns the most nucleophilic character to the -NH2 group in the condensed rings. In Figures S1-S8, see Underlying Data, mailto:https://doi.org/10.7910/DVNs/7KFPUT, Harvard Dataverse, V1 ., the $f^{+}\left(\overrightarrow{r}\right)$ functions for all ligands in the study can be seen.

Figure 15. Fukui function $f^{+}\left(\overrightarrow{r}\right)$ calculated under the FMO approximation $\left({\left|\mathit{\text{LUMO}}\left(\overrightarrow{r}\right)\right|}^2\right)$ for the ligands A) Valacyclovir and B) Penciclovir.

The isovalue was 0.008 in both cases. The figure was created using GaussView 5.0.