Receptor-Based Pharmacophore Modelling of a series of ligands used as inhibitors of the SARS-CoV-2 virus by complementary theoretical approaches, molecular docking, and reactivity descriptors.

System preparation

The receptor structure for the docking experiment was extricated utilizing the following protocols through the crystal Structure of SARS-CoV-2 RNA-dependent RNA polymerase PBD code: 6m71, which was adjusted utilizing the protein preparation wizard module of the Schrödinger suite 2017-1. i) The optimization of the hydrogen bond (H-bond) network and refinement of the protein structure was performed. ii) PropKa utility was utilized to determine the protonation states at physiological pH. iii) The restrained molecular minimization with heavy atoms constrained to a low root-mean-square deviation (RMSD) from the initial coordinates was carried out using the Impact Refinement (Impref) module.^16–18

Alternatively, the molecular structures of the compounds were sketched using Maestro Editor (Maestro, version 11.1, Schrödinger, LLC). Then 3D conformations were gained with the LigPrep module, with ionization/tautomeric states predicted under physiological pH conditions with Epik. Subsequently, energy minimization was used using the protocol with the Macro model using the OPLS2005 force field.

Molecular docking

The docking investigations were carried out using Glide^19,20 with the default parameters and the Standard Precision (SP) model. The docking grid was generated with the co-crystallized ligand at the center using the default settings. To enable the binding of larger ligands, a scaling factor of 0.8 was applied to the van der Waals radii of the nonpolar protein atoms. The Extra Precision (XP) was also employed for induced fit docking (IFD) to allow the protein to undergo backbone, side-chain, or both movements upon ligand docking and to expand alternate receptor conformations suitable for binding ligands with unusual orientations. Finally, considering the extent of residue movement induced by the IFD computation, to evaluate the ligands, we examined the conformations of the most and least active compounds for each molecule in the molecular set. The best predictions of the poses were predicted by 10ns molecular dymane calculations, in order to analyse its stabilization in the active site.

Chemical reactivity analysis

Several previous investigations have established a correlation between quantum similarity and chemical reactivity descriptors.^21–31 Quantum similarity and DFT utilize the density function as an object of study for similarity indexes. Specifically, the Coulomb index can be linked to electronic factors related to chemical reactivity. The global reactivity indices, such as chemical potential (μ),³² hardness (ɳ),³³ and electrophilicity (ω),^34,35 will be calculated using Frontier Molecular Orbitals (FMO) and the energy gap. These chemical reactivity indices (Equations 1-5) provide information about the stability of the systems. Chemical potential measures the tendency of electrons to leave the equilibrium system,³⁶ while chemical hardness measures the resistance of a chemical species to change its electronic configuration.²⁹

The mathematical definition of the electrophilicity index (ω) is related to the stabilization energy of a system when it becomes saturated by electrons from the external environment^34,35:

For this study, the local reactivity descriptors utilized were the Fukui functions. Equations (4, 5) depict how the chemical potential of a system responds to changes in the external field. This is defined as the derivative of the electronic density with respect to the number of electrons at a constant external field.

The nomenclature $f^{+}\left(\overrightarrow{r}\right)$ and $f^{-}\left(\overrightarrow{r}\right)$ have been used to indicate nucleophilic and electrophilic attacks, respectively.^36–38 This methodology employs global and local reactivity descriptors to analyze quantum similarity in the molecular set. All calculations were carried out using the method B3LYP³⁹ and the basis set 6-311xxG(d,p)⁴⁰ which is the result of adding a correction to the 6-311G(d) basis set leading to calculations of electronegativity, hardness, reactivity indices and frontier molecular orbitals and is comparable in quality to those obtained with much larger basis sets (such as Aug-cc-pVQZ and Aug-cc-pV5Z). This method/basis set has been used in combination with the Gaussian 16 package.⁴¹

Molecular analysis: pharmacophore development

To obtain a deeper analysis about the features associated to the biological activity the pharmacophore development has been used using the Zidovudine, Tromantadine, Pyramidine, Oseltamivir, Hydroxychoroquine, Cobicistat, Doravirine (Pifeltro), Dolutegravir, Boceprevir, Indinavir, Truvada, Trizivir, Trifluridine, Sofosbuvir and Zalcitabine. Figure 1 shows the pharmacophore development using the ligands: A) Truvada, B) Trizivir, C) Trifluridine and D) Sofosbuvir. To develop a pharmacophore model, the co-complexed ligand was extracted from its original conformation in the protein and subjected to the pharmacophore development tool available in the Schrödinger suite 2017-1.¹⁶

Figure 1. Pharmacophore development for the ligands analysed.

A) Truvada, B) Trizivir, C) Trifluridine and D) Sofosbuvir.

A receptor-based pharmacophore model was constructed (Figure 1) to identify key features associated with biological activity, namely negative, positive, and aromatic ring. The model incorporated hydrogen bond acceptor, hydrogen bond donor, and hydrophobic features, which showed good agreement with a previously reported model. The validation process yielded an RMSD value of 0.38. Figure 2 shows the molecular docking outcomes for the ligand studied.

Figure 2. Molecular docking outcomes (A) and interaction (B) for zidovudine into de pocket protein crystal structure of SARS-CoV-2 RNA-dependent RNA polymerase PBD code: 6m71.

Figure 2 shows the H-bonds for the best conformation of zidovudine with the residues CYS145 and GLU166. For this ligand the best pose has been with docking score -6.241. The H-bond is formed with the alcohol group with the residue CYS145 and ketone group for the residue GLU166.

As seen in Figure 3, tromantadine showed H-bonds with the residue LEU141 with bond length 1.62Å, GLN189 with bond length 1.54Å and with the residue GLU166 with bond length 1.59Å. The H-bonds with the residues LEU141 and GLN189 are with the amino groups and ketone group for the residue GLU166, like Zidovudine.

Figure 3. Molecular docking outcomes (A) and interaction (B) for tromantadine.

Figure 4 shows the H-Bond with the residues GLY143 with length 1.53Å, GLU166 with length 1.61Å and HIE41 with length 1.71Å. The GLY143 and GLU166 had H-bonds with the ketone group, HIE41 had H-bond with the aromatic ring, and finally GLU166 had a H-bond with an amino group.

Figure 4. Molecular docking outcomes (A) and interaction (B) for cobicistat.

On the other hand, Figure 5 shows the H-bonds for doravirine. This compound had H-bonds with the residues GLY143 with length 1.59Å, CYS145 with length 1.59Å, GLU166 with length 1.47Å and a bond with an aromatic ring with length 1.69Å with the residue HIE41.

Figure 5. Molecular docking outcomes (A) and interaction (B) for doravirine (Pifeltro).

Dolutegravir presented H-bonds with the residues GLU166 with a length of 1.48Å. Other H-bonds were with the residues GLY143 with a length of 1.55Å and GLN189 with a length of 1.73Å (see Table 1).

Finally, Figure 7 shows the docking outcomes for the indinavir ligand. This compound had two interactions with the residue GLU166 with lengths of 1.43Å and 1.53Å. On the other hand, this ligand had a H-bond with the residue HIE41 with a length of 1.68Å.

Electrostatic potential comparison

In the previous section, we have seen that the interactions that form between the ligands and the RNA polymerase can be classified into two main types, purely electrostatic attractions and interactions by delocalisation of charges. In this and the following section, we develop these aspects using DFT calculations. Figure 8 shows the function ESP for the ligands zidavudine, cobicistat and dolutegravir. In the case of zidavudine the ESP on the O2 and O3 atoms has a clear correspondence with the interactions with CYS145 and GLU166 which can be seen in Figure 2. In the case of cobicistat the significant values of the ESP function on the O4 and O7 atoms have a correspondence with the interactions with GLU166 and GLY143 (Figure 4) respectively. Finally, in the case of dolutegravir, the ESP function calculated on the O4 and O7 atoms had a clear agreement with the interactions with CYS145 and GLU166 respectively, which can be seen in Figure 6.

Figure 6. Molecular docking outcomes (A) and interaction (B) for dolutegravir.

Figure 9 shows the Function ESP for the ligands tromantadine, doravirine and indavir. In the case of tromantadine, the ESP on the O1 and N4 atoms have a clear correspondence with the interactions with GLU166 and LEU141 which can be seen in Figure 3. In the case of doravirine the significant values of the ESP function on the O6 and N10 atoms had a correspondence with the interactions with GLU166 and GLY143 (Figure 5) respectively. Finally, in the case of indavir, the ESP function calculated on the N5 atom had a clear correspondence with the interaction with GLU166, which can be seen in Figure 7. In Figures S13-S24 (Extended data⁴³) 10.6084/m9.figshare.22670167.v1, the images of all the ESPs for all the ligands studied can be seen.

Figure 7. Molecular docking outcomes (A) and interaction (B) for indinavir.

Figure 8. Function ESP for the ligands: A) Zidavudine, B) Zidavudine but with the different contributions of each atom to the ESP represented by colours, C) Cobicistat, D) The same as B but for cobicistat, E) Dolutegravir and F) The same as B but for dolutegravir.

The isovalue for A, C and D was -0.04, and the iso. for B, D and E was -0.01. Figures A, C and D were created using GaussView 5.0 and B, D and E using AIMAll (v. 17.11.14).

Figure 9. Function ESP for the ligands: A) Tromantadine B) Tromantadine but with the different contributions of each atom to the ESP represented by colours, C) Doravirine, D) The same as the B but for doravirine, E) Indavir and F) The same as the B but for indavir.

The isovalue for A, C and D was -0.04, and the iso. for B, D and E was -0.01. Figures A, C and D were created using GaussView 5.0 and B, D and E using AIMAll (v. 17.11.14).

Global reactivity descriptors analysis and Fukui function comparison

The study also investigated the global and local chemical reactivity descriptors using DFT calculations. Table 2 presents the calculated global parameters, including chemical potential, chemical hardness, global softness, and global electrophilicity, to compare the chemical reactivity of the ligand sample. As indicated in Table 2, the least reactive molecule is tromantidine, exhibiting the lowest values of electronic chemical potential μ, softness S, and electrophilicity ω (for the chemical potential it refers to its absolute value). On the other hand, it has the highest chemical hardness (η) value. The most reactive compounds are pyramidine and doravirine with the highest values for electronic chemical potential, softness and electrophilicity, as well as the lowest chemical hardness values. The electrophilicity values can have a crucial influence on the stability of the active site of ligands that are stabilized by non-covalent interactions.

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 10 shows the functions for the compounds cobicistat, hydroxychoroquine, indinavir, oseltamivir and tromantadine (A-E respectively), it can be noted that in these five cases the function assigns the most nucleophilic character to a nitrogen atom, mainly to its unshared electron pair. When comparing this figure with Figures 3 and 7 we can see that some important interactions can be justified on this basis; for example, in the case of Indinavir the N5 has an important interaction with GLU166, or in the case of Tromantadine the N4 has an important interaction with LEU141.

Figure 10. Fukui function $f^{-}\left(\overrightarrow{r}\right)$ calculated under the FMO approximation $\left({\left|\mathit{HOMO}\left(\overrightarrow{r}\right)\right|}^2\right)$ for the ligands.

A) Cobicistat, B) Hydroxychoroquine, C) Indinavir, D) Oseltamivir and E) Tromantadine. Isovalue was 0.01 in all cases. The figure was created using GaussView 5.0.

Figure 11 shows that in the case of zidavudine the O3 would have an important interaction as a charge donor with GLU166 (Figure 2). For doravirine, the N10 has an interaction with GLY143 and the N11 with CYS145 (Figure 5). In the case of dolutegravir we have not found any match for the function $f^{-}\left(\overrightarrow{r}\right)$. In Figures S25-S36 (Extended data⁴³), we can see images of all the functions $f^{-}\left(\overrightarrow{r}\right)$ for all the ligands studied.

Figure 11. Fukui function $f^{-}\left(\overrightarrow{r}\right)$ calculated under the FMO approximation $\left({\left|\mathit{HOMO}\left(\overrightarrow{r}\right)\right|}^2\right)$ for the ligands.

A) Zidavudine, B) Doravirine and C) Dolutegravir. The isovalue was 0.01 in all cases. The figure was created using GaussView 5.0.

Figure 12 shows the functions $f^{+}\left(\overrightarrow{r}\right)$ calculated under the FMO approximation $\left({\left|\mathit{\text{LUMO}}\left(\overrightarrow{r}\right)\right|}^2\right)$ for compounds A) Zidavudine, B) Tromantadine, C) Cobicistat, D) Doravirine, E) Dolutegravir and F) Indavir. In the case of tromantadine, the Fukui function $f^{+}\left(\overrightarrow{r}\right)$ on the N3 and N4 atoms justifies interactions by charge attraction towards these atoms with GLN189 and LEU141 respectively. In the case of dolutegravir, the value of the function on O4 indicates a possible interaction by charge delocalisation with GLU143. For the rest of the ligands no matches for the function $f^{+}\left(\overrightarrow{r}\right)$ were found. In Figures S25-S36 in the supporting information (Extended data⁴³) are images of all the functions $f^{+}\left(\overrightarrow{r}\right)$ of all the ligands studied.

Figure 12. Fukui function $f^{+}\left(\overrightarrow{r}\right)$ calculated under the FMO approximation $\left({\left|\mathit{LUMO}\left(\overrightarrow{r}\right)\right|}^2\right)$ for the ligands.

A) Zidavudine, B) Tromantadine, C) Cobicistat, D) Doravirine and E) Dolutegravir. The isovalue was 0.01 in all cases. The figure was created using GaussView 5.0.