Computational studies on potential new anti-Covid-19 agents with a multi-target mode of action

2.1 Molecular electrostatic potential (MEP) mapping studies

Compounds 1–5 were prepared using the Gaussian 09 software package as previously reported (Mohapatra et al., 2021e). Calculations were performed using the B3LYP-D2 (Grimme, 2006)/6-31G** basis set in the Gaussian 16 program suite (Fig. S2) (Marenich et al., 2009). The electrostatic potential distribution for each molecule was measured by using the following equation (Takayanagi et al., 1996); $$\mathrm{\Phi }\left(r\right)={\int }^{\frac{\rho \left(r\right)}{\left|r-r^{\prime }\right|}}dr$$ where, ρ explains both the nuclear and electronic charge of a molecule, with the integration performed over the molecular volume of the same molecule, and r’ is the atomic position relative to same origin.

All molecular electrostatic potentials were colour-graded, with red and blue corresponding to the most negative and most positive potential, respectively, and green and yellow corresponding to a potential halfway between the two extremes.

2.2 Frontier molecular orbitals (FMO) studies

The TD-DFT (time-dependent density functional theory) calculations were reported by using Gaussian 09 platform (Sahu et al., 2021) with B3LYP method and 6–31 + G (d, p) basis set to identify the chemical reactivity of compounds, by calculating E_HOMO and E_LUMO energies and the [E_HOMO - E_LUMO] energy difference (band gap). This also helped to characterise global reactivity descriptors including chemical hardness (η), electrophilicity index (ω), electronic chemical potential (μ), and electronegativity (χ) for all compounds (Sahu et al., 2021).

2.3 Natural bond orbital (NBO) studies

Natural bond orbital (NBO) analysis was conducted in order to explore the bonding aspects related to compounds 1–5, using a set of anti-bond, localized bond, and Rydberg extra valence orbitals (Mohamed et al., 2018; Monika et al., 2021). NBO studies were used to analyse donor–acceptor interactions for all compounds using a second-order perturbation approach (Fock matrix). The strength of the stabilization energy E⁽²⁾ associated with electron delocalisation among donor (i) and acceptor (j) NBO was deduced using the following equation: $$E^{\left(2\right)}=q_i\frac{F^2\left(i,j\right)}{{\epsilon }_{j-{\epsilon }_i}}$$ where, q_i denotes donor orbital occupancy, F(i, j) implies the off diagonal NBO Fock matrix element and ε_j and ε_i are the diagonal elements, High stabilization energy E⁽²⁾ is indicative of a strong interaction between the electrons of the acceptor–donor NBO, suggesting the presence of a higher conjugation system (Sumrra et al., 2018).

2.4 Molecular docking study

The 3D crystal structures of the human ACE2 receptor (PDB ID:1R42), the SARS-CoV-2 RBD S-protein (PDB ID:6MOJ Chain B), Papain-Like protein (PLpro) (PDB ID:7CJD Chain A) and reverse transcriptase (RdRp) (PDB ID:6YYT Chain A) were obtained from Protein Data Bank (https://www.rcsb.org/pdb) and used as biological targets for the molecular docking study (Xu et al., 2020; Xu et al., 2004; Gao et al., 2021; Hillen et al., 2020). All proteins were prepared by removing water molecules, ions, and heteroatoms, with grid maps created around the respective binding sites using AutoDock Tools and Auto Grid programs. Prior to the docking, the 2D and 3D structures of all ligands were drawn using HyperChem v.8.0 and optimized geometrically using the MM + force field (Hocquet and Langgård, 1998). The ligands and the protein active site residues were allowed to be flexible during the docking process, which was conducted using AutoDock Vina (Trott and Olson, 2010). The active site was chosen around 5 Å of the ligan-protein interaction region. The docking grid dimension was taken to be 25 × 25 × 25 Å and the default grid point is taken to be 0.375 Å. The best pose (with lowest docking score) was used for further studies. All protein–ligand interactions were visualized using Biovia Discovery Studio Visualizer v.16.1.0. Ligand efficiencies (LE) were calculated for each compound (Hopkins et al., 2014) by using the following equation; $$\mathit{LE}=\left(1.37/HA\right)\times pI{C}_{50}=\left(1.37/HA\right)\times pKd.$$ where, HA is the number of heavy (non-hydrogen) atoms, pKd is the negative logarithm to the base 10 of the dissociation constant and pIC₅₀ is the negative logarithm to the base 10 of the half-maximal inhibitory concentration.

2.5 Molecular dynamics (MD) simulation studies

With the help of YASARA molecular modeling v.18.4.24, the stability of the protein–ligand (P-L) complexes was determined (Krieger and Vriend, 2014). The P-L complexes were put in a cubic simulation box (density 0.997 g L⁻¹, TIP3P water model, NaCl 0.9 as counter ions) in aqueous medium, keeping at least a 5 Å buffer space around each complex. The systems were energy-minimized using the AMBER14 force field with a steepest gradient approach (100 cycles). The MD simulation was carried out for 40 ns with a frame capture every 25 ps and trajectories were generated from the data captured.

2.6 Pharmacokinetic, drug-likeness and toxicity prediction studies

The pharmacokinetic properties of compounds 1–5 were predicted using the online admetSAR v1/v2 and SwissADME software. The Lipinski’s properties were tabulated by using physicochemical data obtained from PubChem database. Drug-likeness filters (Ghose, Veber, Muegge and Egan) were also determined to enhance the predictions (Mohapatra et al., 2021b). Toxicity of the compounds was predicted by using ProTox II online server (Banerjee et al., 2018).

3.1 MEP mapping studies

The molecular electrostatic potential maps obtained for compound 1–5 are shown in Fig. 1. All compounds displayed reactive centres with the most negative potential predicted at oxygen atoms. Compounds 1–3 showed a more negative electrostatic potential (with darker red oxygen atoms) compared to compounds 4–5. MEP mapping studies are used to predict the chemical reactivity of molecules, with high electron density regions (prone to electrophilic attack) and low electron density regions (prone to nucleophilic attack) showing a negative and a positive electrostatic potential, respectively (Yadav et al., 2021).

Close

Fig. 1

Molecular electrostatic potential (MEP) maps obtained for compounds 1–5 (a: denotes transparent display and b denotes solid display).

Export to PPT

3.2 FMO studies

The electronic transitions taking place from the HOMO to LUMO orbitals of compounds 1–5 are shown in Fig. 2. Compound 1 displayed the highest E_HOMO-E_LUMO gap (5.349 eV) in comparison with 2 and 5 which had the lowest gaps. As a lower transition gap indicates a higher reactivity/lower stability (Choudhary et al., 2013), the energy difference values obtained suggest that compounds 2 and 5 show the greatest reactivity amongst all compounds. The values obtained for the electronic chemical potential (μ), electronegativity (χ), chemical hardness (η), and electrophilicity index (ω) are presented in Table 1. Compound 1 showed the highest value for chemical hardness and electrophilicity, while compound 5 showed the highest value for chemical potential and electronegativity. These results confirmed that compound 1 was the most stable amongst all compounds as high chemical hardness and low electrophilicity are directly linked to high stability (Choudhary et al., 2013).

Close

Fig. 2

[E_HOMO − E_LUMO] energy band gaps obtained for compounds 1–5 (value in electron volt).

Export to PPT

3.3 NBO studies

The NBO plots corresponding to the carbonyl (C⚌O) bond of compound 1–5 are shown in Fig. 3. The examination of the NBO plot obtained for compound 1 revealed that p_z orbital of both the carbon and the oxygen atoms were involved in σ-bond formation with the carbon and oxygen atom showing a 34.95 % and 65.05 % orbital contribution, respectively. The p_y orbital is further involved in the formation of a π-bond with the carbon and the oxygen atom showing a 32.71 % and 67.29 % orbital contribution, respectively. The orbital contributions obtained for all compounds are listed in Table 2, while the corresponding natural electronic configurations of all carbonyl bonds are listed in Table 3.

Close

Fig. 3

Natural Bond Orbital (NBO) plots obtained for compounds 1–5.

Export to PPT

The total Lewis and non-Lewis contributions for compound 1 and all other compounds (2–5) are illustrated in Fig. 3 and Fig. S3, respectively. Compound 1 exhibited a Lewis type structure (99.931%) with a non-Lewis (Valence-2.12% and Rydberg-0.16%) contribution. Compound 2 had a Lewis type structure (98.905%) with a non-Lewis (Valence-0.931% and Rydberg-0.164%) contribution. Compound 3 exhibited a Lewis type structure (98.928%) with a non-Lewis (Valence-0.917% and Rydberg-0.155%) contribution. Compound 4 had a Lewis type structure (98.959%) with a non-Lewis (Valence-0.888% and Rydberg-0.154%) contribution. Compound 5 exhibited a Lewis type structure (98.144%) with a non-Lewis (Valence-1.6215% and Rydberg-0.2337%) contribution.

The electron donor–acceptor orbitals and interaction stabilization energy values obtained for all compounds are illustrated in Table 4. Strong intramolecular interactions were predicted to occur in all compounds. Compound 1–5 were stabilized predominantly by n to σ* transitions. This included n(O17) → σ*(C12-C16), n(O17) → σ*(C13-C14) (10.80 kcal/mol each) for 1; n (O17) → σ*(C12-C13), n(O17) → σ*(C13-C14) (10.46 and 11.94 kcal/mol, respectively) for 2; n(O15) → σ*(C4-C13), n(O15) → σ*(C13-C14) (11.66 and 10.46 kcal/mol, respectively) for 3; n(O14) → σ*(C1-C2), n(O14) → σ*(C2-C3) (10.26 and 11.35 kcal/mol, respectively) for 4; n(O1) → σ*(C2-O13), n(O13) → σ*(O1-C2), n(O13) → σ*(C2-C3) (22.26, 16.04, 9.08 kcal/mol, respectively) for 5. The latter showed additional π(C–C) and antibonding π* interactions, including π(C26-C28) → π*(C25-C27), π(C26-C28) → π(C29-C30), πC29-C30) → π*(C26-C28), πC25-C27) → π*(C29-C30) (12.07, 9.99, 10.03, 11.3 kcal/mol, respectively). This indicated that it showed better conjugation and greater stability compared to other compounds. These π → π* and n → π*interactions have remarkable impact on crystal packing and charge transfer interactions (Sumrra et al., 2018). The bending angles corresponding to C⚌O bonds (i.e. angle of deviation from the direction of the line connecting the two nuclei centres) for all compounds are reported in Table S1.

3.4 Molecular docking studies

The docking scores (binding energy, BE), ligand efficiency (LE), and Ki values obtained for the best conformational positions of 1–5 with ACE2, RBD, PLpro, RdRp, and detailed intermolecular interactions of each ligand with the target proteins are presented in Tables S2–S9, and depicted in Fig. 4. Compounds 4 and 5 displayed the best results against PLpro, ACE2, and RdRp, while compounds 1, 2, 5 displayed best results against RBD. The BE for compound 5 against ACE2, RBD, PLpro, RdRp found to be −7.549, −7.392, −8.031 and −7.883 kcal/mole respectively which showed better interaction and stability of the ligand against the targets. The Lys441 and Gln442 of ACE2; Gly339 and Ser371 of RBD; Tyr273, Asp164 and Pro248 of PLpro; and Gly327 and Pro328 of RdRp forms H-bond with the ligand at the active site. The results obtained for the binding free energy (ΔGbind) of compound 5 with all target proteins using MM/PBSA (GBSA) are presented in Table S10.

Close

Fig. 4

Intermolecular interactions of compound 5 with (a) PLpro, (b) RdRp, (c) RBD, (d) ACE2.

Export to PPT

3.5 MD simulation studies

During MD simulation studies, two important trajectories - ligand movement and ligand conformation − were analyzed and expressed in terms of RMSD in Å. All protein–ligand complexes analysed revealed RMSD values ≤3.0 Å (Fig. 5a). The trajectories of RMSD of Cα of the backbone of P-L complexes are shown in Fig. 5b. The RMSD values obtained for ACE2-5 and RBD-5 were <2.0 Å while PLpro-5 had a value >2 Å. The average RMSDs of the ligand movement trajectories of ACE2-5, RBD-5, PLpro-5 and RdRp-5 were 3.753, 6.038, 4.009 and 6.186, respectively. These values (>3 Å) suggest a re-orientation of the ligands in each protein cavity. Structural investigation shows the ligands are well within the cavity of the active site. The radius of gyration for all complexes remained within <1 Å of their average values, which suggests that the P-L complexes stayed compact throughout the simulation (Fig. 6). The trajectories of radius of gyration of the P-L complexes of all the ligands suggest compact conformations, reflecting the stability of the complexes (Mohapatra et al., 2021f).

Close

Fig. 5

(a) RMSD in Å of the ligand movement of ACE2-5 (blue, avg. 3.753), RBD-5 (red, avg. 6.038), PLpro-5 (green, avg. 4.009), and RdRp-5 (brown, avg. 6.186) complexes during the MD simulations. (b) RMSD in Å of protein–ligand Cα trajectory of ACE2-5 (blue, avg. 1.728), RBD-5 (red, avg.1.319), PLpro-5 (green, avg. 2.715), and RdRp-5 (brown, avg. 1.809) complexes during the MD simulations.

Export to PPT

Close

Fig. 6

Radius of gyration (Rg) in Å of ACE2-5 (blue, avg. 25.236), RBD-5 (red, avg. 18.435), PLpro-5 (green, avg. 18.435), and RdRp-5 (brown, 29.74) complexes during the MD simulations.

Export to PPT

3.6 Pharmacokinetic, drug-likeness and toxicity parameters

The pharmacokinetic and drug-likeness parameters for the ligands (Fig. S4) were calculated and are reported in Tables S11–S16. All investigated compounds obeyed the rules for oral drug-likeness, including the Lipinski’s rule of five (Lipinski, 2000; Lipinski, 2004). In the case of compounds 1–5, the molar refractivity was <100. The topological polar surface area (TPSA) for these compounds was < 70 Å², indicating their potentiality as favourable drug molecules (Dariya and Nagaraju, 2020). Their number of rotatable bonds (≤2) suggested flexibility. These compounds were soluble and highly absorbable in gastrointestinal (GI) tract. Their synthetic accessibility value was nearly 5, indicating the feasibility of their synthesis (Mohapatra et al., 2021b). The toxicity parameters including LD50 and toxicity class of the five compounds are listed in Table S13. Compounds 2 and 4 with lower LD50 values of 750 and 1000 mg/kg predicted to be more toxic. Compounds 1, 3 and 5 having LD50 values of 5000 mg/kg, 4600 mg/kg and 2000 mg/kg are less toxic and can be considered for oral consumption.