Targeting virus-host interaction by novel pyrimidine derivative: an in silico approach towards discovery of potential drug against COVID-19.

The entire human population over the globe is currently facing appalling conditions due to the spread of infection from coronavirus disease-2019 (COVID-19). The spike glycoprotein of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) present on the surface of the virion mediates the virus entry into the host cells and therefore is targeted by several scientific groups as a novel drug target site. The spike glycoprotein binds to the human angiotensin-converting enzyme-2 (hACE2) cell surface receptor abundantly expressed in lung tissues, and this binding phenomenon is a primary determinant of cell tropism and pathogenesis. The binding and internalization of the virus is the primary and most crucial step in the process of infection, and therefore the molecules targeting the inhibition of this process certainly hold a significant therapeutic value. Thus, we systematically applied the computational techniques to identify the plausible inhibitor from a chosen set of well characterized diaryl pyrimidine analogues which may disrupt interfacial interaction of spike glycoprotein (S) at the surface of hACE2. Using molecular docking, molecular dynamics (MD) simulation and binding free energy calculation, we have identified AP-NP (2-(2-amino-5-(naphthalen-2-yl)pyrimidin-4-yl)phenol), AP-3-OMe-Ph (2-(2-amino-5-(3-methoxyphenyl)pyrimidin-4-yl)phenol) and AP-4-Me-Ph (2-(2-amino-5-(p-tolyl) pyrimidin-4-yl)phenol) from a group of diaryl pyrimidine derivatives which appears to bind at the interface of the hACE2-S complex with low binding free energy. Thus, pyrimidine derivative AP-NP may be explored as an effective inhibitor for hACE2-S complex. Furthermore, in vitro and in vivo studies will strengthen the use of these inhibitors as suitable drug candidates against SARS-COV-2. Communicated by Ramaswamy H. Sarma.

Introduction

The world is currently going through a debilitating phase of acute health disaster attributed to the global pandemic brought about by the novel coronavirus disease-2019 (COVID-19) (Enayatkhani et al., 2020; Joshi et al., 2020; Kirchdoerfer et al., 2016; Muralidharan et al., 2020). Sequencing and simultaneous phylogenetic identification of the virus responsible for COVID-19 confirmed that it as a novel β-coronavirus that shared 88% sequence identity with two bat-derived SARS-like coronaviruses (Lu et al., 2020;; Pant et al., 2020; Wang et al., 2020). Additionally, it was shown that this coronavirus (CoV), termed as 2019-nCoV (Elfiky, 2020; Gorbalenya et al., 2020), shared 79.5% sequence identity with SARS-CoV (Elfiky & Azzam, 2020; Lu et al., 2020;; Wang et al., 2020; Xia et al., 2020) which caused the severe acute respiratory syndrome pandemic in 2012. Therefore, this newly identified virus was called as SARS-CoV-2 (Aanouz et al., 2020). The virus gains entry into the host cells via the transmembrane spike (S) glycoprotein that forms a homotrimeric structure extending outwards from the viral surface (Belouzard et al., 2012; White et al., 2008). The monomeric S protein is comprised of two distinct functional domains S1 and S2. The S1 domain is responsible for binding to the host cell receptor, while S2 domain helps in the fusion of the viral and cellular membranes (Jitendra Subhash et al., 2020; Kirchdoerfer et al., 2016). The S protein is cleaved at the boundary between the S1 and S2 subunits thereby adopting a prefusion conformation that remains bound non-covalently to the surface (Walls et al., 2020). The S1 protein subunit contains a unique receptor-binding domain (RBD) in its C-terminal region, which contributes towards the stabilization of the prefusion state of the membrane-anchored S2 subunit containing the fusion machinery (Shang et al., 2020; Walls et al., 2020). Transmembrane protease serine 2 (TMPRSS-2) further cleaves the S protein at a specific site located immediately upstream of the fusion peptide (Elmezayen et al., 2020; Shen et al., 2017). This site specific cleavage induces a set of irreversible conformational changes that essentially activates the protein for membrane fusion (Shen et al., 2017). After binding of the RBD to the hACE2 receptor on target cells, the heptad repeat 1 (HR1) and 2 (HR2) domains in the S2 subunit of S protein interact with each other to form a six-helix bundle (6-HB) fusion core, bringing viral and cellular membranes into close proximity for fusion and infection (Du et al., 2009; Li et al., 2019; Xia et al., 2020; Zhu et al., 2013).

Scientists throughout the world are currently looking for effective and immediate therapeutic strategies to intervene and prevent the infectivity or transmission of SARS-CoV-2. Till date, no suitable vaccines or anti-viral agents have been put forward.

Several types of pyrimidine bases like thymine, cytosine and uracil are found in nature, which forms the framework or building blocks of genetic materials, i.e. DNA and RNA (Sharma et al., 2014; Sohrabi et al., 2018). This prime importance makes pyrimidine moiety suitable for broad therapeutic applications (Shakibapour et al., 2019; Sharma et al., 2014). Pyrimidine moieties have gained considerable attention in the chemotherapy of various diseases like cancer (Xie et al., 2011), HIV, diabetes, cardiovascular, bacterial, fungal diseases (Sharma et al., 2014), and several viral diseases (Balzarini & McGuigan, 2002). It is also used as anti-inflammatory drugs (Amir et al., 2007), analgesics (Vega et al., 1990), antipyretic drugs, anti-leishmanial drugs, herbicidal agents and anti-oxidants (Abu-Hashem et al., 2011). As pyrimidine base is a constituent of both DNA and RNA; hence, it is useful in the chemotherapy of both DNA and RNA viruses (Sharma et al., 2014). Compounds with N-heterocyclic scaffolds are highly effective against a diverse range of diseases, and it has already been studied for their potential pharmacological activities. In particular, pyrimidine derivatives have demonstrated outstanding biological activities. For example, Dayvigo (Lemborexant), Inrebic (Fedratinib) are a few examples of drugs containing pyrimidine scaffolds which received FDA-approval in 2019. Because of the enormous therapeutic importance, we selected some well-characterized pyrimidine substituted phenols (Kumar & Rao, 2018; Table S1, ) to investigate their efficiency in binding to the interface of the hACE2-S protein complex and modulate the pattern of infectivity of the SARS-CoV-2 virus. In this study, we aim to identify diaryl pyrimidine derivatives as a potential lead molecule which may bind at the interface of the hACE2-S protein complex with high affinity. Our computational analyses show that AP-NP, AP-3-OMe-Ph and AP-4-Me-Ph have the potential to bind at the interface of the hACE2-S protein complex with high affinity. Thus, the further experimental studies with these lead molecules may facilitate the anti-COVID drug development processes.

Table 1.

The binding energy and the amino acids involved within specific distance of diaryl pyrimidine derivatives against hACE2-S protein complex.

Compound	Binding affinity (kcal/mol)	Interacting amino acid residue of hACE2
AP-NP	−8.95	LEU 391, PHE 390, LEU 73, TRP 69, PHE 40, ASP 350, ASN 394, ARG 393
AP-Ph-4-I	−7.9	PHE 40, ASP 350, TRP 349, ALA 348, ASP 382, TYR 385, PHE 390
AP-3-OMe-Ph	−8.1	PHE 40, TRP 349, ALA 348, ASP 382, TYR 385, ARG 393, PHE 390, ASP 350, GLY 352
AP-4-Me-Ph	−8.1	ALA 348, TRP 349, SER 47, SER 44, PHE 40, ASP 350, TYR 385, HIS 401, ASP 382
AP-Ph-4-OMe	−7.8	TYR 385, ASP 382, HIS 401, ALA 348, TRP 349, PHE 40, PHE 390, GLY 352, ARG 393
AP-Ph-4-Br	−7.9	PHE 40, TRP 349, ALA 348, ASP 350, ASP 382, TYR 385, PHE 390
Chloroquine	−5.7	TRP 69, PHE 40, ASN 394, ARG 393, PHE 390, LEU 391, ALA 99, LEU 100, LEU 73

Materials and methods

Preparation of ligands and receptor

The three-dimensional structure of the hACE2-S protein complex (human angiotensin-converting enzyme 2 and spike glycoprotein complex, PDB ID: 6VW1) was downloaded from the RCSB protein data bank (Shang et al., 2020). The structure of all pyrimidine derivatives was prepared in Discovery studio 2020. The conformation of all ligands was refined using ‘Clean geometry’ command and the structure with minimum clean energy was saved in pdb format.

Molecular docking of diaryl pyrimidine derivatives on hACE2- S protein complex

The molecular docking analysis of the AP-NP compound was performed on the X-ray crystal structure of the hACE2-S protein complex (PDB ID: 6VW1). First, polar hydrogen and gasteiger charges were added on the spike glycoprotein and AP-NP compound using MGL tool (Morris et al., 2009) and saved in pdbqt format. The molecular docking tool AutoDock Vina (Trott & Olson, 2010) was used for the docking of AP-NP on the hACE2-S protein complex. To determine the probable binding site, blind docking was performed. For this, the hACE2-S protein complex was kept rigid, while AP-NP was kept flexible to allow exploration of probable binding sites. Both molecules were covered by grid box with a dimension of 116 Å × 82 Å × 86 Å with grid spacing 1 Å. Five sets of docking were performed with exhaustiveness of 100. The AutoDock Vina produced 9 docked conformations with binding energies in the range of −7.95 to −8.95 kcal/mol. Among these conformations, the conformation, where AP-NP was observed to bind at the proximal helix site present at the junction of the interacting hACE2 and S1 domains, possessed the lowest binding free energy. Therefore, the helical region was chosen for local docking. For this, the helical region was covered by the grid box with a dimension of 58 Å × 82 Å × 120 Å with grid spacing 1 Å and docking was done with exhaustiveness of 100. Five sets of local docking were performed, and the conformation with the lowest binding energy was chosen for analysis. The analyses of the results/outputs from AutoDock Vina were done using MGL tools 1.5.6 (Morris et al., 2009), and intermolecular interactions were determined using PyMoL (DeLano, 2002). The molecular docking experiments of the other diaryl pyrimidine derivatives were performed in a similar fashion as done with the AP-NP.

Molecular Dynamics (MD) Simulation study

Using GROMACS v5.1.4 biomolecular simulation package, with CHARMM27 force field and TIP3P water model, all-atoms molecular dynamics (MD) simulations were carried out on the coordinates of hACE2-S protein complex docked with compounds: AP-NP, AP-3-OMe-Ph and AP-4-Me-Ph (Abraham et al., 2015; Mokaberi et al., 2019; Sharifi-Rad et al., 2020). The topology and parameter files for ligands parameters were repaired as define by Zoete et al. (2011), Mongre et al. (2019), and Sharifi-Rad et al. (2020). The simulation box was prepared with 10 Å buffer distance, and the protein–ligand complex was placed at the centre of box, padding around with water molecules. Along with counter ions for neutralizing the system, the salt concentration of 0.15 M (0.15 M NaCl) was used for the MD simulation (Joung & Cheatham, 2008). PBC condition was defined for all three x, y and z directions (Darden et al., 1993). All energy minimization and simulations were performed at physiological temperature, 300 K (Kumar et al., 2019; Mishra et al., 2018). Energy minimization was done estimating the steepest descent, followed by the conjugant gradients involving 50,000 steps for each. SHAKE algorithm was applied to constrain all the bonds involving hydrogen atom. Particle Mesh Ewald (PME) was applied for long-range electrostatic forces. The equilibration of system was done in two steps, NVT and NPT for the period of 500 ps. Pressure and temperature during the simulation were maintained by Berendsen thermostat (Berendsen et al., 1987) and Parrinello-Rahman pressure (Parrinello & Rahman, 1980), respectively. LINC algorithm was applied to constrain the bonds and angles (Hess et al., 1997). LJ potential with a cut off of 0.10 nm was defined for Van der Waals interactions. Finally, the production run was performed on NPT ensemble for a period of 100 ns, and the trajectory was updated at the time interval of 10 ps to record the energy and velocity. All production runs were performed on CUDA enabled Tesla GPU machine (DELL T640 with V100 GPU), and OS Centos 7 (Prakash et al., 2018; Singh et al., 2019) and the GROMACS utilities were used for the analyses of obtained MD trajectories.

Binding free energy estimation

The binding free energy of the protein–ligand complexes was calculated using MM-PBSA (Molecular Mechanics-Poisson-Boltzmann Surface Area), which describes the structural stability, spatial orientation and molecular interactions of ligands in the active site of the protein (Batt et al., 2012; Prakash & Luthra, 2012; Sastry et al., 2013). The binding free energy abbreviated as ΔGbinding can be defined as, where Gcomplex represents the total free energy of protein-ligand complex, Gprotein as the free energy of protein, Gligand used as the free energy of ligand and < > represents the ensemble average.

Neglecting the entropy term (TS), the binding free energy can be approximately written as where ΔEMM is the change in the molecular mechanics interaction energy (gas-phase) upon ligand binding calculated as the sum of the changes in the bonded energy, electrostatics and Van der Waals interactions upon ligand binding. ΔGsolv is the change in solvation free energy upon ligand binding.

Furthermore, ΔGsolv can be written as where is the change in the polar part of solvation free energy and is the change in nonpolar part of solvation free energy as a result of ligand binding to the proteins. In this work, Poisson–Boltzmann (PB) method was used for the estimation of the polar part of the solvation free energy.

In the current work, binding free energy () was calculated utilizing the MMPBSA.py script of the AMBER package (Wang et al., 2016). Single trajectory protocol approach was utilized for the computation, which involves only the simulation of complex forms. An ionic strength of 0.15 M was used for the PBSA calculations. Considering the convergence issues associated with the MM-PBSA calculation, only last 20 ns data was used.

ADME analysis for drug likeliness

The drug likeliness property of all diaryl pyrimidine derivatives was monitored as described earlier (Jayaram et al., 2012; Lipinski, 2004). The pdb structure of all diaryl pyrimidine derivate compounds was constructed, and ADME properties, as explained in Lipinski’s rule of five, i.e. molecular weight, solubility, H-bond donor, H-bond acceptor and molar refractivity were calculated by online software tool.

Results and discussion

Molecular docking analysis

Designing the targeted drug molecules by using cheminformatics tools is the initial step in drug discovery (Xu & Hagler, 2002). In this study, our focus was to investigate the efficacy of diaryl pyrimidine derivative compounds in obstructing the interaction between spike glycoprotein of COVID-19 and the hACE2 receptor, which is considered to be an essential step towards the progress of infection. The interaction of spike protein and hACE2 facilitates the attachment and internalization of the virus in the host cell (Belouzard et al., 2012). The COVID-19 virus attaches itself to the host cell using S1 domain spike glycoprotein, whereas the S2 domain of spike protein facilitates the internalization of the virus into the host membrane to begin infection (Belouzard et al., 2012; Du et al., 2009). So both domains play a crucial role in the initialization of the infection in the host cell (Enmozhi et al., 2020). The study presented here, using molecular docking tool, revealed significant binding interaction of diaryl pyrimidines with the interface of the hACE2-S receptor complex (hACE2-spike protein complex: PDB ID 6VW1) (Table 1).

In this study, we have analyzed the binding energies of all diaryl pyrimidine derivatives on the hACE2-S complex. Here we found that all diaryl pyrimidine derivatives bind at the interface of the hACE2-S protein complex. Interestingly, these molecules share the same binding site, which lies near the binding junction of hACE2 and C-terminal S1 domain of spike glycoprotein. Among these diaryl pyrimidine derivatives, AP-NP, AP-3-OMe-Ph and AP-4-Me-Ph displayed high binding affinity towards the hACE2 receptor (Figure 1) with the binding energies of −8.95, −8.1 and −8.1 kcal/mol, respectively (Table 1). These pyrimidine containing phenols binds with the active site through hydrogen bonding using –OH group of phenol part and an –NH2 group of the amino pyrimidine ring. In addition to that, the non-covalent interactions such as π–π interaction are also involved in the stabilization of the ligand within the active site.

Figure 1.

Binding site of AP-NP, AP-3-OMe-Ph and AP-4-Me-Ph on hACE2-S protein complex. (a) The putative binding site of AP-NP, AP-4-Me-Ph and AP-3-OMe-Ph on hACE2-S protein complex. The C-terminal S1 domain of spike protein is shown in orange colour, hACE2 is shown in dark green colour. (b) The zoomed view of binding site of AP-NP on hACE2-S protein complex. (c) The interacting amino acid residues of hACE2 with AP-NP. The AP-NP is shown in cream colour. (d) The zoomed view of binding site of AP-4-Me-Ph on hACE2-S protein complex. (e) The interacting amino acid residues of hACE2 with AP-4-Me-Ph. The AP-4-Me-Ph is shown in purple colour. (f) The zoomed view of binding site of AP-3OMe-Ph on hACE2-S protein complex. (g) The interacting amino acid residues of hACE2 with AP-3-OMe-Ph. The AP-3-OMe-Ph is shown in light pink colour.

In particular, the significant binding energy observed in case of AP-NP can be attributed to such non-covalent hydrophobic interactions of phenyl as well as napthyl ring with non-polar amino acid residues such as LEU 391, PHE 390, LEU 73, TRP 69 and PHE 40 along with the interaction through hydrogen bonding with polar amino acid residues such as ASP 350, ASN 394 and ARG 393 (Figure 1(b, c)). In case of AP-3-OMe-Ph, PHE 40, PHE 390, TRP 349 of hACE2 forms hydrophobic interaction with both the aryl rings attached with pyrimidine scaffold (Figure 1(d, e)). In addition to that, hydrogen bonds were also observed with ALA 348, ASP 382, ARG 393, ASP 350, and GLY 352. In case of AP-4-Me-Ph, ALA 348, TRP 349, SER 47, SER 44, PHE 40 of hACE2 forms hydrophobic interaction with the toluene ring. However, the phenol-pyrimidine base forms hydrogen bonding with ASP 350, TYR 385, HIS 401 and ASP 382 (Figure 1(f, g)).

The similarity in the binding energy of AP-3-OMe-Ph with AP-4-Me-Ph could be explained due to the presence of tetrahedral –OMe group instead of the planar aromatic group such as naphthyl. It also has the same hydrophobic interaction as AP-4-Me-Ph. This might be due to their structural similarity. Compared to the AP-4-Me-Ph, additional hydrogen bonding interactions were observed with ARG 393 residue due to the involvement of oxygen atom of –OMe group.

CQ, a well-known drug that is suggested for use in the COVID-19 therapy, was reported to interfere with the terminal glycosylation of spike protein of SARS-CoV-2S and thereby reducing the initial infection (Hu et al., 2020; Vincent et al., 2005). Hence, we also used CQ as a positive control and performed molecular docking experiment. Similar to all diaryl pyrimidine derivatives, we found that CQ also binds at the interface of the hACE2-S protein complex (). We observed that CQ has lesser binding affinity (−5.7 kcal/mol) with hACE2-S complex, which may be due to its high affinity towards its primary target, i.e. endosome (Al-Bari, 2017).

These selected inhibitors generated through the molecular docking studies were projected onto further analyses for molecular simulation dynamics and estimation of the free energy of binding using MM-PBSA.

To examine the spatial stability and mechanistic aspects of conformational dynamics underlying the molecular interaction of diaryl pyrimidine derivatives, AP-NP, AP-3-OMe-Ph and AP-4-Me-Ph with hACE2-S protein complex, we performed MD simulation in an aqueous environment for the period of 100 ns, at physiological temperature (300 K). We analyzed the MD trajectory by employing structural parameters, RMSD (root mean square deviation) and Rg (radius of gyration), H-bond (hydrogen bonding), the minimum distance between protein and ligands, principal component analysis (PCA) and free-energy landscape (FEL) which revealed the dynamic stability of ligand-bound complex (Kumar et al., 2019; Mishra et al., 2018; Sarma et al., 2020). MM-PBSA was used for the estimation of binding free energy, which describes the conformational stability of ligand at the binding site, molecular interactions and the binding affinity with protein.

To determine the global conformational stability of ligands docked with hACE2-S protein, we computed all atoms Cα-backbone RMSD with respect to the initial structure, as shown in Figure 2. The result shows that the RMSD plot of hACE2-S protein complex docked with AP-NP quickly achieves the stable equilibrium during the initial few nanoseconds and stable conformational dynamics is observed till the end of the simulation.

Figure 2.

The evolution plot of Cα RMSD in water at 300 K of hACE2-S protein docked complex with compounds, AP-NP (black), AP-3-OMe-Ph (red) and AP-4-Me-Ph (green).

The RMSD trajectory of AP-3-OMe-Ph shows a continuous rise in RMSD during the initial 0–25 ns. The initial rise in RMSD suggested the spatial arrangement of the active site to accommodate the ligand. The ligand adjusted well at the binding groove; i.e. the protein–ligand complex achieved equilibrium at ∼25 ns and it remained stable for the period of 100 ns. Similarly, we observed consecutive two drifts in the RMSD trajectory of hACE2-S protein docked complex with AP-4-Me-Ph, during the initial 0–20 ns. The system attains equilibrium only after 20 ns, and a stable conformation of the protein–ligand complex propagated for the remaining simulation time. However, the consequences of small drifts during ∼40–70 ns indicate relatively lower conformational stability of AP-4-Me-Ph as compared to AP-NP and AP-3-OMe-Ph. Another structural order parameter, Rg, is commonly used to determine the structural compactness and stability of protein during the simulation. Thus, the time evolution plot of Rg is calculated, as shown in Figure 3.

Figure 3.

Time evolution plot of radius of gyration (Rg) in water at 300 K of hACE2-S protein docked complex with compounds, AP-NP (black), AP-3-OMe-Ph (red) and AP-4-Me-Ph (green).

In this figure, we can see that the conformation dynamics of the hACE2-S protein-AP-NP complex remain compact and stable with an average Rg value of 2.7 ± 0.01 nm, throughout the simulation period. We also did not observe any significant structural changes in the Rg trajectory of AP-3-OMe-Ph, and the protein–ligand complex remains stable with Rg value 4.4 ± 0.01 nm. However, the Rg plot of AP-4-Me-Ph shows two small drifts of ∼0.02 nm. Although, it is well settled around 35 ns and the progression of stable complex structure for the remaining period of simulation. Thus, RMSD and Rg suggest that all three compounds remain bound to the active site during the simulation period with minimal perturbation; however, AP-NP appeared more stable complex as compared to AP-3-OMe-Ph and AP-4-Me-Ph.

To further confirm the molecular interaction stability of ligands with protein, we also measured the occupancy of H-bond interactions between protein-ligands during the simulation. Besides maintaining the conformational architecture and stability of protein structure, H-bonds play a critical role in the stability of the protein–ligand complex. Molecular docking results showed the presence of two H-bonds between AP-NP and hACE2-S protein, which is consistently observed during the simulation time (Figure 4). However, during 0–20 ns and 60–100 ns, only one H-bond remains consistent, and the appearance of two H-bond is seen only for the period of 20–60 ns.

Figure 4.

Time evolution plot of hydrogen bonds (H-bond) between hACE2-S protein and ligands, AP-NP, AP-3-OMe-Ph and AP-4-Me-Ph. The representation of plots color as shown in Figure 1.

The molecular interaction of AP-3-OMe-Ph shows the occupancy of maximum three H-bonds, but only two remain stable throughout the simulation period. Appearance of three H-bond seen only between 40 and 60 ns. The maximum occupancy of four-five H-bonds is seen with AP-4-Me-Ph as observed in molecular docking; however, only two remains consistent during the simulation. But, considering the occurrence of appearance and disappearance of H-bond interactions with ligands, substantially only two H-bonds with AP-NP and AP-3-OMe-Ph observed persistent for the maximum period of simulation.

To maintain the dynamic stability of ligand at the binding pocket of protein, consistency of the spatial orientation of ligand and the optimal distance from the binding cavity of protein plays a significant role in the interaction stability of the protein-ligand complex. Thus, the time evolution plot of the mean distance of all three ligands monitored during the simulation (Figure 5).

Figure 5.

Time evolution plot of minimum (min) distances from the center of binding grooves and ligands. Color scheme as shown in Figure 1.

From Figure 5, it can be clearly depicted that the mean distance of all three ligands remains in the range of 3.7–4.2 Å, which is observed consistent almost throughout simulation time. Only, AP-4-Me-Ph moved out the cavity approximately at ∼15 ns and ∼70 ns, however, it appeared transiently, and all three ligands continue as well occupied to the binding pocket of hACE2- S protein. Furthermore, the lower flexibility of average distance with AP-NP again provides the elegance evidence of stable molecular interaction with the target protein.

Essential dynamics

Essential dynamics (ED) govern the coherent changes in the conformational dynamics of protein–ligand complex while the simulation, which can be calculated by projecting the first two principal components (PCs), PC1, and PC2 along the native structure. By diagonalizing the covariance matrix of eigenvectors, we computed PCs which determine the overall combined motion of the Cα-atoms, and localization of PCs in conformational space correlated to the protein function. ED plot of all three complexes is shown in Figure 6.

Figure 6.

Essential dynamics of hACE2-S protein with AP-NP (black), AP-3-OMe-Ph (red) and AP-4-Me-Ph (green). Conformational ensembles sampled with the projection of principal components PC1 and PC2.

From Figure 6, it is observed that the collective motion of AP-NP confined to small conformation space, reflecting stable collective atomic fluctuations. Differently, AP-4-Me-Ph explores broader conformational space during the simulation, depicting the larger collective motion of atoms; thus, it indicates a relatively less stable structure. We can see that the collective motion of AP-3-OMe-Ph shows broader conformation space than AP-NP; however, it is significantly restricted to less occupied space as compared to AP-4-Me-Ph. Thus, ED analysis indicates a better affinity of AP-NP towards ACE2-S-protein than AP-3-OMe-Ph and AP-4-Me-Ph, respectively.

Estimation of the FEL

We also performed FEL analysis to examine the minimum energy conformational ensembles of all three protein–ligand complexes, using the first two PCs (PC1 and PC2) of the protein backbone atoms as reaction coordinates (Kumar et al., 2018; Sarma et al., 2020). This allows an estimation of the free energy difference of conformational stability of the hACE2-S protein bound complex with AP-NP, AP-3-OMe-Ph, and AP-4-Me-Ph with high accuracy (Figure 7).

Figure 7.

Free energy landscape calculated using the principal components PC1 and PC2 as reaction coordinates of hACE2-S protein with AP-NP, AP-3-OMe-Ph and AP-4-Me-Ph.

From Figure 7, we can see that AP-NP displays less conformational mobility, restricting to a more confined conformational space within single energy minima. This result displayed that molecular interactions of AP-NP contributed to the stabilization of the ligand-bound complex of hACE2-S protein; thus, the complex structure remains connected through the low energy barriers. The broader and sallow energy basin of AP-3-OMe-Ph characterized that ligand-bound complex navigated through the wide range of conformations; however, with small excursions, conformational ensembles shifted the dynamic equilibrium which stayed stable and confined to the low energy basin. Whereas, the sallow and fragmented energy basin of ligand AP-4-Me-Ph suggests that the conformational ensemble of complex structure broadly experiences the two different conformational spaces. The complex structure readily moves out the native basin, highly flexible structure navigated the wide conformational space, before entering to the stable energy basin. This result signifies that the proper positioning of AP-4-Me-Ph orientation in the binding grooves may require conformational adjustment, hence, unstable distributions of the binding configuration ensemble of complex transverse to stable energy basin through the higher energy barrier. Thus, the FEL plot suggests that the compound AP-NP almost stayed constant at the bonding grooves of hACE2-S protein and remains well stabilized through the molecular interactions.

Estimation of the binding free energy

The binding free energy of the protein–ligand complex provides a robust approximation of molecular interactions, which includes bonded and non-bonded interactions and the effect of solvent underlying the protein–ligand binding affinity (Kumar et al., 2019; Pandey et al., 2017). Bonded and non-bonded terms are defined by the van der Waals and electrostatic interactions. Whereas, the solvation free energies in terms of polar and non-polar solvation free energy. This analysis was performed on the fully converged trajectory of protein–ligand for the last 20 ns data (Pandey et al., 2018; Wang et al., 2016). The computed binding free energy (ΔGbinding) and the contribution of other energy parameters on the molecular interaction stabilities of AP-NP, AP-3-OMe-Ph and AP-4-Me-Ph are summarized in Table 1. Results show that the lowest binding energy, ΔG = −30.89 ± 2.24 kcal mol−1, was obtained for ligand AP-NP, suggesting the highest binding affinity with hACE2-S protein. Out of the three compounds, AP-4-Me-Ph possess the highest binding free energy −25.29 ± 2.14 kcal mol−1, which indicates the lowest binding affinity towards the target protein complex. Whereas, the binding free energy of AP-3-OMe-Ph (−28.39 ± 2.69 kcal mol−1) was observed slightly higher than AP-NP, thus, considered as the ligand with a moderate range of binding affinity. Furthermore, the results show a larger contribution of the Van der Waals energies as compared to others like electrostatic energies suggesting that ligands at the binding pocket are predominantly stabilized by hydrophobic interactions. Thus, our analyses provided substantial evidence that hydrophobic amino acid residues at the groove of the hACE2-S-protein complex is crucial for the structure-based drug development in the therapy of COVID-19.

Lipinski’s rule of five

To determine whether any compound with a particular biological activity has the potential to serve as a pharmacological agent/drug, the Lipinski rule of five is generally used. This rule acts as a filter in in-silico analyses and helps to screen potential drugs during the very beginning of the drug designing process, thus minimizing the cost of exercises, time, and labour in case of clinical drug development (Gombar et al., 2003; Hughes et al., 2011). The rule judges some of the basic molecular properties a compound possesses, like absorption, distribution, metabolism and excretion (ADME) for a selected range (Hughes et al., 2011). If a compound possesses at least two of the properties, it’s good to go (Jayaram et al., 2012; Lipinski, 2004). The potential diaryl pyrimidine derivatives used in this study were found to pass all the five criteria mentioned in Lipinski’s rule (Table 2). So it also supported the therapeutic importance of the selected diaryl pyrimidine moieties.

Table 2.

ADME values for selected diaryl pyrimidine derivatives derivatives.

S. No.	Compound/ Ligand	ADME properties (Lipinki’s rule of five)		Drug likeliness
		Properties	Values
1.	AP-NP	Molecular weight (<500 Da)	313	Yes
		LogP (<5)	4.3
		H-bond donor (5)	3
		H-bond acceptor (<10)	4
		Molar refractivity (40–130)	96.5
		Violations	NO
2.	AP-3-OMe-Ph	Molecular weight (<500 Da)	293	Yes
		LogP (<5)	3.1
		H-bond donor (5)	3
		H-bond acceptor (<10)	5
		Molar refractivity (40–130)	85.5
		Violations	NO
3.	AP-4-Me-Ph	Molecular weight (<500 Da)	277	Yes
		LogP (<5)	3.4
		H-bond donor (5)	3
		H-bond acceptor (<10)	4
		Molar refractivity (40–130)	83.7
		Violations	NO

Conclusion

In-silico studies provide the best alternatives to screen potential drug candidates. In this study, we have investigated the ability of some diaryl pyrimidine derivative compounds for binding to hACE2-S-protein complex and obstructing its receptor binding ability using molecular docking. Our docking results were supported by other analyses estimating the free energy of binding using MM-PBSA, free energy landscape, and dynamics simulation. We found that all the chosen compounds were capable of binding to the interface of human cell receptor ACE2 and spike protein and with high affinity and fulfill all the criteria of Lipinski’s rule of five. We believe that all of these compounds can be effective anti-COVID drugs and, therefore, should be experimentally verified by researchers working on this scheme.

Click here for additional data file.