Identification of some novel oxazine substituted 9-anilinoacridines as SARS-CoV-2 inhibitors for COVID-19 by molecular docking, free energy calculation and molecular dynamics studies.

Coronavirus disease (COVID-19), a life-threatening disease, is caused by SARS-CoV-2. The targeted therapeutics of small molecules helps the scientific community to fight against SARS-CoV-2. In this article, some oxazine substituted 9-anilinoacridines (A1-A48) was designed by docking, MM-GBSA and molecular dynamics (MD) simulation studies for their COVID-19 inhibitory activity. The docking of ligands A1-A48 against SARS-CoV-2 (PDB ID: 5R82) are performed by using Glide module, in silico ADMET screening by QikProp module, binding energy using Prime MM-GB/SA module, MD simulation by Desmond module and atomic charges were derived by Jaguar module of Schrodinger suit 2019-4. Compound A38 has the highest G-score (-7.83) when compared to all the standard compounds which are proposed for COVID-19 treatment such as ritonavir (-7.48), lopinavir (-6.94), nelfinavir (-5.93), hydroxychloroquine (-5.47) and mataquine (-5.37). Compounds A13, A23, A18, A7, A48, A46, A32, A20, A1 and A47 are significantly active against SARS-CoV-2 main protease when compared with hydroxychloroquine and mataquine. The residues GLN19, THR24, THR25, THR26, LEU27, HIE41, SER46, MET49, ASN119, ASN142, HIE164, MET165, ASP187, ARG188 and GLN189 of SARS-CoV-2 main protease play a crucial role in binding with ligands. The in silico ADMET properties of the molecules are within the recommended values. The binding free energy was calculated using PRIME MM-GB/SA studies. From the ligands A38, A13, A23, A18, A7, A48 and A46 with significant Glide scores may produce significant COVID-19 activity for further development. Compound A38 was subjected to MD simulation at 100 ns to study the dynamic behaviour of protein-ligand complex.Communicated by Ramaswamy H. Sarma.

Introduction

Recently, coronavirus disease 2019 (COVID-19), a life-threatening disease, spreads throughout the world. The coronaviruses (CoVs) belong to severe acute respiratory syndrome (SARS) and the Middle East respiratory syndrome (MERS) viruses (De Wit et al., 2016; Song et al., 2019), which contain single positive-stranded RNA. It was affected first in China and quickly spread in most of the countries (Holshue et al., 2020; Wang et al., 2020; Zhou, Hou, et al., 2020). According to the WHO data, as on third week of May 2020, more than 5 million people in the world are affected by COVID-19; out of these, more than 3.3 lakh people are died. Many of the cases affected by COVID 19 are found to be asymptomatic (Gu et al., 2020; Holshue et al., 2020; Lu et al., 2020; To et al., 2020), it is worthy of consideration, the detail current evidence and understanding of the transmission of SARS-CoV, MERS-CoV and SARS-CoV-2 and discuss pathogen inactivation methods on CoVs is very important (Chang et al., 2020; Huang et al., 2020; Huang & Herrmann, 2020; Zhang et al., 2020; Zhou, Yang, et al., 2020).

The CoV transfers from human to human and spreads more vigorously. In this pathetic situation, it is important to discover some novel drugs to treat COVID-19. As COVID-19-positive cases are increased day by day more aggressively, the finding of novel drugs for COVID-19 treatments is one of the most encouraging tasks for target-based treatment.

As a part of our ongoing research on searching the potent biological molecules against various diseases by in silico and wet laboratory methods (Kalirajan et al., 2011; Kalirajan, Kulshrestha, et al., 2012; Kalirajan, Mohammed Rafick, et al., 2012; Kalirajan et al., 2013; Kalirajan, Gowramma, et al., 2017; Kalirajan, Sankar, et al., 2017; Kalirajan, Kulshrestha, et al., 2018; Kalirajan, Mohammed Rafick, et al., 2018; Kalirajan, Gaurav, et al., 2019; Kalirajan, Pandiselvi, et al., 2019), we have designed and evaluated various heterocyclic compounds for their biological activities. Using different modules (Glide, QikProp, Prime and Desmond) of Schrödinger suite LLC, various computational methods such as molecular docking, ADMET screening, binding free energy calculations and molecular dynamics (MD) simulations were performed to find the interactions responsible for COVID-19 inhibition. These studies will provide the requirement of key structural features in the design of potential drug candidates.

It is evident that from the literature studies, 9-anilinoacridines have various pharmacological activities such as anticancer (Kapuriya et al., 2008; ; Wakelin et al., 2003), antimicrobial (Nadaraj et al., 2009), antioxidant (Dickens et al., 2002), antimalarial (Anderson et al., 2006), analgesic (Sondhi et al., 2002), antileishmanial (Giorgio et al., 2007), antinociceptive (Llama et al., 1989), acetyl cholinesterase inhibitors (Recanatini et al., 2000), antiherpes (Goodell et al., 2006) and so forth. Amsacrine that is 9-anilinoacridine derivative was a primary DNA-intercalating agent which is more significant and the chromophore intercalates with DNA base pairs (Rastogi et al., 2002). Similarly, oxazine derivatives were also reported for various biological activities (Kalirajan et al., 2009; Kalirajan, Muralidharan, et al., 2012) such as antimicrobial, anti-cancer, and larvicidal. As a part of our ongoing research work, upon looking for various biological activities, we have designed 9-anilinoacridine analogues bearing the oxazine build-up by docking studies with Schrodinger suit-2019-4. The recently designed 9-anilinoacridines showed significant hindrance with COVID-19.

Materials and methods

Molecular docking studies

The 3D crystal structure of COVID-19 protein called SARS-CoV-2 main protease receptor co-crystallized with 6-(ethylamino) pyridine-3-carbonitrile (PDB ID: 5R82, Resolution: 1.31 Å) was retrieved from the protein data bank. The protein was prepared using protein preparation wizard of epic module (Sastry et al., 2013) of Schrödinger suite 2019. The protein structure is a monomer, which was prepared by removing similar binding sites, unnecessary water molecules and also refining bond orders. Missing chain atoms are added by using the prime module of Schrödinger suite 2019. Protein minimization was performed using optimized potentials for liquid simulations-3 (OPLS-3) molecular force field with root mean square difference (RMSD) of crystallographic heavy atoms kept at 0.3 Å. A grid box was generated to define the centroid of the active site.

The designed ligands (A1–A48) were prepared by using LigPrep module of Schrodinger suite 2019-4. 2D structures were converted to 3D structures, as well as energy minimization and optimized for their geometry, desalted and corrected for their chirality. The ionization and tautomeric states were generated between pH 6.8 and 7.2 by using Epik module. The ligands A1–A48 were minimized using OPLS-3 force field in Schrodinger suit 2019-4 until a RMSD of 2.0 Å was achieved. A single low energy ring confirmation per ligand was generated, and the optimized ligands were used for docking analysis.

The compounds are prepared by LigPrep module of Schrödinger suite and docked in to a catalytic pocket of COVID-19 by using Glide module of Schrödinger suite 2019 in XP (extra precision) mode (Jacobson et al., 2004). The binding modes with best glide G scores were selected. These scores perceive positive lipophilic, hydrogen-bonding and metal–ligand associations and punish steric conflicts. The XP visualizer of Glide module was used to analyze the results of docking studies (Friesner et al., 2006). The Glide score of the molecules was compared with the Glide score of the standard compounds ritonavir, lopinavir, nelfinavir, hydroxychloroquine and mataquine (Barros et al., 2020; Cao et al., 2020) which is recommended for COVID-19. The Gide scoring capacity is mainly dependent on docking parameters like lipophilic perseverance in which the ligands are covered in the lipophilic pocket. The electrostatic forces and hydrogen bonding with ligands are other parameters to increase the binding affinity. Some negative parameters like XP penalties, rotational penalties etc. will diminish the Glide score.

The in silico ADMET properties of the proposed ligands A1–A48 were determined by QikProp module of Schrodinger suit. The prepared ligands were subjected to the QikProp module to determine the ADMET properties such as molecular weight, dipole moment, number of H-bond donors and acceptors, log P values, violence in rule of five, and oral absorption.

The atomic charges were derived using B3LYP hybrid functional combined with 6-31G basis set in Jaguar module of Schrodinger (Govindarajan et al., 2012).

Binding free energy calculation by using prime/MM-GBSA approach

The binding free energy of ligand-receptor complex and post docking energy minimization studies were performed using prime molecular mechanics-generalized born surface area (MM-GB/SA) of Schrödinger 2019. The energy for minimized XP docked pose of ligand receptor complex was calculated using the OPLS3 force field and generalized-born/surface area (GB/SA) continuum VSGB 2.0 solvent model (Li et al., 2011).

MD simulation

The stability of the docked A38/5R82 complex was investigated by performing a 100 ns MD (Guo et al., 2010) simulation study. The complex in the explicit solvent system with OPLS3 force field was studied using Desmond module of Schrödinger 2019-4. The molecular system was solvated with crystallographic water (TIP3P) molecules (Jorgensen et al., 1983) under orthorhombic periodic boundary conditions for 10 Å buffer region. The overlapping water molecules are deleted, and the system was neutralized by adding Na+ as counter ions. The total system consists of 39,270 atoms and 11,507 water molecules. An ensemble (NPT) of Nose-Hoover thermostat (Martyna et al., 1992, 1994) and barostat was applied to maintain the constant temperature (300 K) and pressure (1 bar) of the systems, respectively. A hybrid energy minimization algorithm with 1000 steps of steepest descent followed by conjugate gradient algorithms was utilized. Another algorithm, limited memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) algorithm with convergence threshold gradient of 1 kcal/mol/Å was also employed for energy minimization. A Smooth Particle Mesh Ewald method for calculating long range electrostatic interactions with a cut-off radius of 9 Å for short range van der Waals and Coulomb interactions was used. Multiple time step RESPA integration (reference system propagator algorithms) was used in the dynamics study for bonded, near and far-bonded interactions with 2, 2 and 6 fs, respectively. The data were collected for every 100 ps, and the obtained trajectory was analyzed with Maestro graphical interphase.

Results and discussion

The results are summarized in Tables 1–3 and Figures 1–6. The results revealed that the COVID-19 inhibitory property of the compounds A1–A48 greatly depended on the chemical nature of the substituents. The chemical structures of oxazine substituted 9-anilinoacridines are given in Figure 1.

Table 1.

Docking studies for oxazine substituted 9-anilinoacridines A1–A48 with COVID-19 (5R82).

Cpd	Glide score	Lipophilic EvdW	H Bond	XP Electro	Low MW	XP Penalties	Rot Penal
A38	−7.829	−5.174	−1.081	−1.98	0	0	0.146
A13	−6.78	−5.77	−1.48	−0.57	0	0	0.16
A23	−6.65	−6.14	−0.61	−0.29	0	0	0.14
A18	−6.37	−5.75	−1	−0.36	0	0	0.15
A7	−6.12	−5.83	−1.04	−0.41	0	0	0.16
A48	−5.991	−5.129	−0.882	−0.195	−0.138	0	0.202
A46	−5.978	−5.145	−1.007	−0.248	−0.232	0	0.172
A32	−5.805	−5.499	−1.218	−0.289	0	0	0.164
A20	−5.80	−5.27	−1.28	−0.35	0	0	0.12
A1	−5.79	−5.77	−0.1	−0.11	0	0	0.13
A47	−5.784	−5.096	−0.454	−0.239	−0.185	0	0.161
A36	−5.478	−5.628	−0.649	−0.398	−0.005	1	0.171
A5	−5.34	−4.85	−0.89	−0.31	0	0	0.17
A12	−5.33	−5.02	−0.7	−0.5	0	0	0.15
A19	−5.22	−5.24	−0.53	−0.04	0	0	0.15
A25	−5.15	−5.33	−0.64	−0.29	0	0	0.14
A11	−5.05	−5.36	−0.7	−0.32	0	0	0.16
A45	−4.99	−5.261	−0.53	−0.176	−0.022	0	0.174
A8	−4.82	−5.42	−0.35	−0.05	0	0	0.16
A35	−4.80	−5.316	−0.171	−0.132	0	0	0.147
A9	−4.72	−5.46	−0.21	−0.05	0	0	0.16
A22	−4.56	−5.38	−0.04	−0.06	0	0	0.12
A40	−4.50	−5.217	−0.036	−0.093	0	0	0.134
A30	−4.48	−5.42	−0.082	−0.065	0	0	0.153
A31	−4.44	−5.683	−0.83	−0.292	0	1	0.153
A2	−4.40	−5.38	−0.08	−0.04	0	0	0.13
A16	−4.38	−5.24	−0.08	0.01	0	0	0.18
A17	−4.34	−5.14	−0.19	−0.06	0	0	0.13
A34	−4.321	−5.231	0	−0.039	0	0	0.147
A4	−4.23	−5.22	−0.12	−0.05	0	0	0.17
A21	−4.22	−5.06	0	−0.08	0	0	0.12
A15	−4.22	−5.32	0	0.01	0	0	0.18
A29	−4.202	−5.202	0	−0.063	0	0	0.153
A6	−4.16	−5.33	−0.04	−0.01	0	0	0.17
A24	−4.04	−5.31	−0.55	−0.31	0	1	0.13
A37	−3.96	−5.211	−0.223	−0.285	−0.058	1	0.182
A10	−3.89	−4.92	0	0	0	0	0.16
A3	−3.86	−4.56	−0.7	−0.31	0	0	0.13
A14	−3.82	−3.94	−0.7	−0.45	0	0	0.16
A39	−3.80	−5.169	−0.024	−0.065	0	0	0.179
A26	−3.07	−2.29	−1.35	−0.16	0	0	0.12
A44	−3.04	−5.171	−0.072	−0.017	−0.022	0	0.174
A33	−2.98	−4.977	0	−0.068	0	0	0.164
A27	−2.92	−4.943	0	0.044	−0.025	0	0.175
A43	−2.86	−5.088	−1.048	0.039	−0.022	1	0.174
A28	−2.67	−4.414	0	−0.133	0	0	0.155
A41	−2.58	−3.308	−0.7	−0.278	0	1	0.139
A42	−2.23	−3.137	0	−0.14	0	0	0.139
Hydroxychloroquine (Std)	−5.47	−3.15	−1.75	−0.69	−0.38	0.5	0
Ritonavir (Std)	−7.48	−6.57	−1.7	−0.46	0	0	0.22
Lopinavir (Std)	−6.94	−6.01	−1.12	−0.26	0	0	0.24
Nelfinavir (Std)	−5.93	−4.86	−1.18	−0.49	0	0	0.24
Metaquine (Std)	−5.37	−5.3	−0.7	−0.38	−0.06	0	0.18
6-(ethylamino)pyridine-3-carbonitrile (Bound inhibitor)	−4.40	−2.9	−0.7	−0.19	−0.5	0	0.29

Figure 1.

Structures of compounds A1–A48.

The docking studies of the ligands to protein active sites were performed by an advanced molecular docking program Glide module of Schrodinger suite 2019 Maestro-12.2 version for determining the binding affinities of the compounds. The designed analogues were docked towards the COVID-19 (PDB ID: 5R82) in order to ascertain their inhibitory activity. The analogues show a best fit RMSD value of 0.2. As shown in Table 1, it is clearly demonstrated that Compound A38 has the highest G-score (−7.83) when compared to all the standard compounds which are proposed for COVID-19 treatment such as ritonavir (−7.48), lopinavir (−6.94), nelfinavir (−5.93), hydroxychloroquine (−5.47) and mataquine (−5.37). Compounds A38, A13, A23, A18, A7, A48, A46, A32, A20, A1, A47 and A36 are significantly active against COVID-19 with Glide score more than −5.5 when compared to currently used drug hydroxychloroquine (−5.47) and mataquine (−5.37). The above compounds have good affinity to the receptor due to more lipophilic character and hydrogen bonding.

The docking results of the compounds exhibited similar mode of interactions with COVID-19 and the binding pocket of the residues between GLN19 and GLN189. The results are summarized in Table 1. Almost all Compounds A1–A48 with COVID-19 (PDB ID: 5R82) are docked in the same binding pocket. From Figure 2a, the amino acid residues GLN19, THR24, THR25, THR26, SER46, HIE41, ASN142, HIE164 and GLN189 are making polar region and the amino acids LEU27, MET49 and MET165 are making hydrophobic interaction with the ligand. The hydroxyl group of the phenyl ring is interacted by hydrogen bonding with ASN 119. The Glide scores are mainly increased due to the lipophilic evidence of the aromatic moieties. The 2D-ligand interaction diagram of compounds A38, A13, A23, A18, A7, A48, A46 and A32 with COVID-19 (PDB ID: 5R82) are given in Figure 2(a–h).

Figure 2.

(a) Ligand Interaction of compound A38 with COVID-19 (5R82). (b). Ligand Interaction of compound A13 with COVID-19 (5R82). (c) Ligand Interaction of compound A23 with COVID-19 (5R82). (d) Ligand Interaction of compound A18 with COVID-19 (5R82). (e) Ligand Interaction of compound A7 with COVID-19 (5R82). (f) Ligand Interaction of compound A48 with COVID-19 (5R82). (g) Ligand Interaction of compound A46 with COVID-19 (5R82). (h) Ligand Interaction of compound A32 with COVID-19 (5R82).

From the obtained binding modes, it was illustrated that the ligands formed hydrophobic, hydrogen bonding and other interactions with different residues GLN19 to GLN189 surrounding the active pocket. Ligand A13 exhibited hydrogen bonding interaction with ASN142 (H-bond length: 3.06 Å) and ASN119 (H-bond length: 1.79 Å) residues and these are shown in Figure 3. The presence of aromatic features and different heterocyclic rings majorly contributed towards lipophilic factors (Figure 4).

Figure 3.

Hydrogen bonding interaction of compound A13 with COVID-19 (5R82).

Figure 4.

Hydrophilic/lipophilic map of compound A23 with COVID-19 (5R82).

The ADMET properties for compounds A1–A48 can be investigated by the in silico method by using QikProp module of Schrödinger suite 2019-4. The molecular weight of the compounds is between 380 and 544 g mol−1. The dipole moment of the compounds is between 1 and 5.7. The estimated number of hydrogen bond donors of the compounds is in the range of 3–4. The estimated number of hydrogen bond acceptors of the compounds is in the range of 4.5–6.75. The number of likely metabolites of the compounds is in the range of 2–6. The number of violations of Lipinski’s rule of five is 0–2. The percentage of human oral absorption for the compounds is between 70% and 100%; Therefore, almost all the ADMET properties of the compounds are within the recommended values. The results of the ADMET properties for compounds A1–A48 are shown in Table S1 (supplementary material).

The atomic charges were derived using B3LYP hybrid functional combined with 6-31G basis set in Jaguar module of Schrodinger. The HOMO, LUMO and molecular surface property (ESP mean) values are illustrated in Table S3 (supplementary material).

Molecular docking was additionally assessed with MM-GBSA free restricting vitality, which is identified with the post scoring approach for COVID-19 (PDB ID: 5R82) target and the values are shown in Supplementary Table S2 (supplementary material). Prime MMGBSA DG bind, the binding free energy, is calculated with the equation:

From the results of MM-GB/SA studies, the ΔG bind values were observed in the range of −26.01 (cpdA35) to −63.24 kcal/mol (cpdA14) for significantly active compounds and also dGvdw values, dG lipophilic values and the energies are positively contributing towards total binding energy. The stability of docking complex will be predicted by MMGBSA scoring function. The Glide score and MM-GBSA free energy are obtained by the docking of ligands into the coupling pocket, which are more stable.

A 100-ns molecular dynamic simulation was performed to understand the molecular insights involved in binding of A38 in the active pocket of 5R82. From the obtained trajectory analysis, the RMSD of protein Cα, backbone and heavy atoms was observed in the range of 1.25–2.1 Å and 1.45–2.45 Å, respectively (Figure 5). The Cα atoms fluctuated in the range of 1.25–2.1 Å and finally stabilized after 80 ns of simulation with a RMSD value of 1.5 Å. Higher fluctuations (up to 2.1 Å) in RMSD were observed at 39–42 ns and 79 ns. Over the course of simulation, stable hydrophobic interactions were observed with His41 and Met165. The flexibility of residues on ligand binding is analyzed using metrics of root mean square fluctuations (RMSFs) (Figure 6). High flexibility was observed with RMSF of Ser1 (Cα: 4.67 Å; backbone: 4.77 Å) and Val303 (Cα: 5.95 Å; backbone: 6.02 Å), whereas low RMSF was reported for residues Leu205 (Cα: 0.37 Å; backbone: 0.38 Å) and Thr175 (Cα: 0.43 Å; backbone: 0.46 Å). The stability of ligand with respect to protein was indicated using ligand RMSD where a moderate fluctuation of 1.72 Å was observed (Figure 6).

Figure 5.

RMSD (Å) of simulated protein 5R82.pdb in complex with inhibitor A38 during 100 ns MD simulation.

Figure 6.

RMSF of simulated protein 5R82.pdb in complex with inhibitor A38 during 100 ns MD simulation.

A total of 27 ligand contacts (Figure 7) were formed with amino acids of protein, from Gln19 to Met49, Asn119, Phe140 to Pro168 and Phe181 to Gln192. From Figure 8, it was stated that the ligand is stabilized by forming majority of hydrophobic interactions (5%–90% of simulation time) with residues His41, Met49, Cys145, Met165, Leu167, Pro168 and Phe181. It also formed hydrogen bonds with Gly19, Asn28, His41, Asn119, Asn142, Gly143, Ser144 and Gln189 over the course of 10%–65% of simulation trajectory. The timeline representation of protein-ligand contacts is outlined in Figure 8.

Figure 7.

Interaction fraction of 5R82.pdb complexed with A38 during 100 ns MD simulation.

Figure 8.

Time line representation showing different contacts formed by A38 in complex with 5R82.pdb during 100 ns MD simulation.

The 2D-trajectory interaction diagram (Figure 9) depicts that the hydrogen bond formed by the docking pose with Asn119 is preserved in the MD trajectory pose. The methoxy group in compound A38 donated one hydrogen bond to Asn119 with 24%, the nitrogen of acridine ring accepted one hydrogen bond with Arg188 and Ser144 through water molecule with 21% and 51% of total simulation time, respectively. The nitrogen of the oxazine ring accepted one hydrogen bond from Gly143 and the amino group accepted hydrogen bonds from Gly143 and Ser144 with 58% and 51%, respectively, of total simulation time. The pi-pi stacking of acridine with His41 was 67%. The RMSF (Figure 10) of ligand with respect to initial frame was observed in the range of 0.5–1.48 Å. Other ligand properties (Figure 10) such as radius of gyration, molecular surface area, solvent accessible surface area and polar surface area of ligands were observed in the range of 5.85–6.07 Å, 444.94–453.42 Å2, 228.75–453.338.82 Å2 and 149.94–164.73 Å2, respectively.

Figure 9.

2D interaction diagram of A38 in complex with 5R82.pdb during 100 ns MD simulation.

Figure 10.

Ligand properties of A38 in complex with 5R82.pdb during 100 ns MD simulation.

Conclusion

In conclusion, from the oxazine-substituted 9-anilinoacridines A1–A48, many compounds have significant binding affinity with SARS CoV-2 main protease. Molecular docking and binding free energy calculation studies were performed to find the possible binding modes of ligands and the influence of favorable and non-favourable interactions within the active pocket of COVID-19 protein. MD simulation for the highly active inhibitor A38 in complex with protein 5R82 revealed that the stabilization of ligand was achieved due to the formation of uninterrupted hydrophobic interactions. The results demonstrated for further modifications in pharmacophoric features may help in improvement of inhibitory activity. The in silico structuring strategy embraced in the present investigation helped for recognizing some lead molecules and furthermore may somewhat clarify their useful impact for further determinations like in vitro and in vivo assessments. The results from the in silico study exhibited that compounds A38, A13, A23, A18 and A7 may be significantly active against COVID-19 with remedial possibilities and are probably going to be helpful after further refinement.

Click here for additional data file.