Analogue discovery of safer alternatives to HCQ and CQ drugs for SAR-CoV-2 by computational design.

COVID-19 outbreak poses a severe health emergency to the global community. Due to availability of limited data, the selection of an effective treatment is a challenge. Hydroxychloroquine (HCQ), a chloroquine (CQ) derivative administered for malaria and autoimmune diseases, has been shown to be effective against both Severe Acute Respiratory Syndrome (SARS-CoV-1) and SARS-CoV-2. Apart from the known adverse effects of these drugs, recently the use of CQ and HCQ as a potential treatment for COVID-19 is under flux globally. In this study, we focused on identifying a more potent analogue of HCQ and CQ against the spike protein of SAR-CoV-2 that can act as an effective antiviral agent for COVID-19 treatment. Systematic pharmacokinetics, drug-likeness, basicity predictions, virtual screening and molecular dynamics analysis (200 ns) were carried out to predict the inhibition potential of the analogous compounds on the spike protein. This work identifies the six potential analogues, out of which two compounds, namely 1-[1-(6-Chloroquinolin-4-yl) piperidin-4-yl]piperidin-3-ol and (1R,2R)-2-N-(7-Chloroquinolin-4-yl)cyclohexane-1,2-diamine interact with the active site of the spike protein similar to HCQ and CQ respectively with augmented safety profile.

Introduction

Coronavirus Disease-2019 (COVID-19), triggered by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), rapidly became a global pandemic and was reported from more than 200 countries. The outbreak of COVID-19 emerges as a serious risk to public health [1]. As of January 12, 2021, there have been more than 89 million cases of COVID-19 globally, including 1,940,352 deaths reported (https://covid19.who.int/). Such a massive number of infections and fatalities, calls for an imperative development of an effective treatment and affordable therapeutics to recover from the pandemic. FDA approved the use of chloroquine (CQ) and hydroxychloroquine (HCQ) as antimalarial and also for autoimmune diseases, for the treatment of COVID-19 patients as Emergency Use Authorization (EUA) with cautions issued afterwards [[2], [3], [4]].

CQ and HCQ are analogous chemical compounds whose mechanisms as a weak base and immunomodulatory properties are the same. Several studies have suggested the role of CQ in inhibiting viral entry, replication, and post-translational modification [[5], [6], [7], [8], [9], [10], [11]]. Two different mechanisms have been reported on the role of CQ in inhibiting viral entry. These drugs act on the biosynthesis of sialic acid by inhibiting quinone reductase 2 [12,13] and simultaneously interacting with angiotensin-converting enzyme 2 receptors (ACE2), and impacting glycosylation [14,15]. CQ also interrupts the replication of the virus by increasing endosomal pH, thereby inhibiting virus-endosome fusion [16]. CQ also interferes with proteolytic processes and hampers post-translational modifications of viral proteins [17]. Similarly, HCQ also interferes with glycosylation, endosomal fusion, and lysosomal activity by modulating pH [18]. HCQ also impacts the signaling of Toll-like receptors (TLRs) by increasing endosomal pH. Immunomodulatory effects of HCQ suppresses major histocompatibility complex (MCH) class II expression by inhibiting T cell activation, expression of CD145, and release of cytokines [19].

An excessive dose of CQ reportedly causes acute poisoning and death [20]. HCQ, an analogue of CQ, developed with a hydroxyl group into CQ was much less toxic than its parent compound [21]. Recent reports have indicated that the administration of HCQ also resulted in potential adverse effects on dermatomyositis patients and caused retinal toxicity [22,23]. Due to the adverse effects of HCQ and CQ and taking into consideration the treatment, there is a necessity to develop a potentially safer analogue of these drugs with no or minimal side effects. We developed an in silico framework for the discovery of analogue compounds and to identify potentially safer alternatives of CQ and HCQ using virtual screening and molecular dynamics simulations. Virtual screening was carried out in this study by the generation of potent analogues and molecular docking of the selected ligand library against the crystal structure of SARS-CoV- 2 spike protein that exhibits inhibitory binding as antiviral, drug-likeness and safety profiling to treat COVID-19 infection and progression.

Computational details

Generation of analogous library

CQ and HCQ have been reported to have a possible inhibitory effect on the entry of SAR-CoV-2 in the host cell. Structural analogues of these two quinolone molecules were collected from the ZINC drug-like database [24] of the webserver SwissSimilarity [25]. To derive a potential alterative of these compounds, we used the SwissSimilarity web-based tool (http://www.swisssimilarity.ch/), which offers the possibility to perform Ligand Based Virtual Screening. The small molecules libraries cover drugs, bioactive compounds, commercially available molecules and virtual libraries. The molecules most similar to the reference compound are searched, along with the similarity score. The similarity-based virtual screening infers that all the chemicals extracted from the database are similar in the structure to the compound under query and have comparable biological properties vis-à-vis activities [26]. Tanimoto coefficient (T) is the most standard similarity measure for matching chemical structures. Screened structures are generally considered similar if their Tanimoto coefficient is greater than eighty-five percent. Their three-dimensional structure was downloaded in the structure data format (SDF) from the PubChem database [27], which was then converted into Protein Data Bank (PDB) file format using the Open Babel (http://openbabel.org/wiki/Main_Page).

Drug-likeness prediction/assessment

The SwissADME, an open-source tool, uses the Simplified molecular-input line-entry system (SMILE) to calculate various pharmacological properties. The drug-likeness screening of the Swiss ADME includes six rule-based methods such as Lipinski rule-of-five, Ghose filter, Veber Rule, Egan rule and Muegge filter which were applied for the selection of potent analogue for CQ and HCQ [28].

ADMET prediction

ADMET (Absorption, Distribution, Metabolism, Excretion & Toxicity) is vital in assessing the pharmacodynamics activities of the analogue. Several computational tools are available to predict and analyze ADMET profiles of the drug-like molecules. In this study, the following pharmacokinetics properties [29] namely, Blood-Brain Barrier (BBB), Human Intestinal Absorption (HIA), P-glycoprotein inhibitor, CYP450 inhibitory promiscuity (CYP) inhibitor, and carcinogenicity of the ligands were predicted using admetSAR server [30] (http://lmmd.ecust.edu.cn/admetsar2/).

Protein target selection

Among all the 29 viral protein targets found in SAR-CoV-2, the spike protein was selected as it is found to have a potential role in the entry of the virus inside the host [31]. The crystal structure of the novel coronavirus spike receptor-binding domain complexed with its receptor Angiotensin-Converting Enzyme 2 (ACE2) (PDB: 6LZG, 2.5 Å) [32] was obtained from the Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB-PDB) [33]. SARS-CoV-2 spike protein is a heterodimer, consisting of chains A and B. For the molecular docking experiment, the hACE2 protein (chain A) was deleted and chain B from the PDB file of the complex was used. Also, all non-standard residues including water moieties were removed using Discovery Studio [34].

Molecular docking

Molecular docking was performed using AutoDockTools-1.5.6 [35]. Polar hydrogens and Kollman charges were added to the protein and a pdbqt format file was saved. The pdbqt files of each ligand were generated after adding all hydrogen and Kollman charges to each structure. The grid box was set for protein as 90 x 90 x 90 with 0.653 Å spacing to perform a blind (full coverage of the spike protein and complexes) docking simulations. Ten Lamarckian Genetic Algorithm runs were carried out for each ligand. For all the other parameters of docking, default settings were used. All the computational runs were carried out on Cygwin [36] for the generation of both grid parameter file (.gpf file) and a docking parameter file (.dpf file) for each specific ligand. The docked conformations of each ligand were selected based on binding affinity and the top-ranking conformations. The best-selected poses were then visualized using BIOVIA Discovery Studio [34].

Molecular dynamics (MD) simulation

MD simulations of the potent lead compounds, CQ1 and HCQ1 along with their parent drugs CQ and HCQ were carried out using GROMACS 5.1.4 [37]. The spike glycoprotein topology file was prepared by using the GROMOS96 43a1 force field. The ligand (CQ, CQ1, HCQ, HCQ1) topology files were built using the PRODRG server [38]. These complexes were solvated using SPC water models in a dodecahedron periodic box with boundaries extending nearly 10 Å in all directions. The counter ions were added to neutralize the total charge of the system. The steepest descent algorithm followed by conjugate gradient protocol for 50,000 steps with a cut-off value of 1000 kJ·mol−1 was applied for energy minimization [39]. Equilibration of the four complexes was performed in two steps: NVT and then NPT for a period of 100ps each. Temperature coupling with a V-rescale was performed at a temperature of 300 K and a time constant of 0.1 ps, and pressure coupling was done with a Berendsen bath with a time constant 2.0 ps [40]. The Linear Constraint Solver (LINCS) algorithm was used to constrain the bond lengths of heavy atoms [41]. Finally, a 200 ns production run was performed. The coordinate trajectories were saved every 2 ps for the entire time-period [42]. The molecular dynamics analysis was conducted using gmx rms, gmx rmsf, gmx gyrate, gmx sasa, gmx covar, and gmx anaeig module of the Gromacs package.

The principal components analysis of the spike-protein was performed by diagonalizing and solving the eigenvalue (represents the magnitude of motion along with the direction) & eigenvectors (represents of the direction of the motion) for the covariance matrix [43] to obtain functional motions for bio-molecules. Free energy landscape (FEL) represents possible conformations taken by a protein in MD simulation and the Gibbs free energy was calculated using probability distribution from the essential plane composed of the first two eigenvectors. The free energy is plotted along the two order parameters such as RMSD and Rgyr [44] to obtain FEL of the complexes simulated in this investigation. gmx sham tool was used for the construction of FEL.

Basicity prediction

Basicity and acidity of the drugs determine their functional role in targeting endosomal acidification in the biological system, particularly in modulating pH of CQ and HCQ for inhibiting viral replication [45]. The geometry optimizations of CQ and HCQ and their respective analogues were performed using density functional theory (DFT) B3LYP with 6-31 + G (d, p) basis set using Gaussian 16 suite [46]. In the current study, global reactivity descriptors such as ionization potential, electron affinity, chemical potential (μ), electronegativity (χ), chemical hardness (η), and electrophilicity (ω) were calculated [47,48]. The basicity and acidity parameters [49,50] were predicted as mentioned in equations (Eq.1), (Eq. 2)) to understand the reactivity of the selected compounds as compared to those of CQ and HCQ. Where η = chemical hardness; χ = electronegativity.

Results and discussion

We collected a total of 800 analogues of CQ and HCQ from the ZINC database of the webserver SwissSimilarity and performed similarity ranking analysis. Based on the Tanimoto coefficient analysis, 32 analogues (28 CQ and 4 HCQ) were selected with a similarity score above 0.85 for further screening. Details of these analogues of CQ and HCQ with their similarity score, ZINC ID, IUAPC name, and chemical structures are listed in Table 1 and respective SMILEs are given in the supplementary data of Table 1.

Table 1

Analogues of Hydrochloroquine and Chloroquine with their similarity score.

Compounds	IUPAC Name & ZINC Id	PubChem CID	Tanimoto Coefficient	Chemical Structures
HCQ1	1-[1-(6-Chloroquinolin-4-yl)piperidin-4-yl]piperidin-3-ol	72876338	0.886	Image 3
	ZINC95367069
HCQ2	1-(7-chloroquinolin-4-yl)piperidin-4-ol	60631468	0.883	Image 4
	ZINC40412048
HCQ3	2-[4-(7-Chloroquinolin-4-yl) morpholin-2-yl] ethanamine	45177282	0.871	Image 5
	ZINC32509033
HCQ4	[1-(7-Chloroquinolin-4-yl)piperidin-3-yl ] methanol	60413266	0.869	Image 6
	ZINC40412383
CQ1	1R,2R)-2-N-(7-Chloroquinolin-4-yl)cyclohexane-1,2-diamine	93453548	0.946	Image 7
	ZINC38050616
CQ2	(1S,2S)-2-N-(7-chloroquinolin-4-yl)cyclohexane-1,2-diamine	57227818	0.95	Image 8
	ZINC38050615
CQ3	N'-(7-chloroquinolin-4-yl)-N-cyclohexylethane-1,2-diamine	224506	0.957	Image 9
	ZINC01683221
CQ4	N-(4-aminobutyl)-7-chloroquinolin-4-amine	11770817	0.946	Image 10
	ZINC19721342
CQ5	(1R,2S)-2-N-(7-chloroquinolin-4-yl)cyclohexane-1,2-diamine	93453549	0.952	Image 11
	ZINC38050617
CQ6	N'-(7-Chloroquinolin-4-yl)-N-ethylpropane-1,3-diamine	273882	0.933	Image 12
	ZINC01709131
CQ7	N-(7-Chloroquinolin-4-yl)-N′,N′-dimethylbutane-1,4-diamine	3728380	0.984	Image 13
	ZINC01729567
CQ8	N4-(7-chloroquinolin-4-yl)-n1-methylpentane-1,4-diamine	225111	0.989	Image 14
	ZINC02042606
CQ9	Didesethyl Chloroquine	122672	0.974	Image 15
	ZINC06020572
CQ10	3-N-(7-chloroquinolin-4-yl)-1-N,1-N-dimethylbutane-1,3-diamine	4065492	0.959	Image 16
	ZINC06176159
CQ11	N'-(7-Chloroquinolin-4-yl)-N-(2-methylpropyl)propane-1,3-diamine	223166	0.946	Image 17
	ZINC01706242
CQ12	Desethyl Chloroquine	95478	0.993	Image 18
	ZINC02042694
CQ13	N'-(7-chloroquinolin-4-yl)-N-cyclohexylpropane-1,3-diamine	225119	0.975	Image 19
	ZINC01596768
CQ14	1-(7-Chloroquinolin-4-yl)-N,N-dimethylpiperidin-3-amine	60291042	0.974	Image 20
	ZINC78617773
CQ15	N-(7-chloroquinolin-4-yl)-N′,N′-diethylpropane-1,3-diamine	3805581	0.962	Image 21
	ZINC01542123
CQ16	(3R)-1-(7-chloroquinolin-4-yl)azepan-3-amine	96578389	0.947	Image 22
	ZINC82133910
CQ17	N-(7-chloroquinolin-4-yl)-N′-propan-2-ylpropane-1,3-diamine	223165	0.945	Image 23
	ZINC01706241
CQ18	1-(7-chloroquinolin-4-yl)azepan-3-amine	70760365	0.945	Image 24
	ZINC82133908
CQ19	N-(7-Chloroquinolin-4-yl)-N′,N′-di(propan-2-yl)ethane-1,2-diamine	11962135	0.919	Image 25
	ZINC20552561
CQ20	N-(7-chloroquinolin-4-yl)-N′,N′-dimethylpropane-1,3-diamine	11608635	0.915	Image 26
	ZINC01542122
CQ21	N'-(7-Chloroquinolin-4-yl)-N,N-diethylethane-1,2-diamine	408190	0.907	Image 27
	ZINC00006792
CQ22	1-N-(7-Chloroquinolin-4-yl)-2-N,2-N-dimethylpropane-1,2-diamine	11507234	0.906	Image 28
	ZINC37985881
CQ23	(1-Quinolin-4-ylpiperidin-4-yl) methanamine	61785796	0.904	Image 29
	ZINC44514898
CQ24	N-(6-Chloroquinolin-4-yl)-N′,N′-diethylpropane-1,3-diamine	408567	0.899	Image 30
	ZINC01596764
CQ25	N-(7-Chloroquinolin-4-yl)-N′-methylpropane-1,3-diamine	6472985	0.88	Image 31
	ZINC06746756
CQ26	(E)-N-(7-Chloroquinolin-4-yl)-N′,N′-diethylbut-2-ene-1,4-diamine	10093137	0.878	Image 3
	ZINC33956413
CQ27	7-Chloro-N-(5-pyrrolidin-1-ylpentan-2-yl)quinolin-4-amine	220624	0.869	Image 33
	ZINC08579986

Drug-likeness involves various structural features and molecular properties that determine how similar a particular ligand molecule is to the known drugs. These filters are important in drug discovery to design derivatives/analogues within the purview of drug-likeness space [51]. From the similar compounds identified, the 4 analogues of HCQ and 25 out of 28 analogues of CQ do not violate any of the drug-likeness rules. Three analogues of CQ (ZINC01596768, ZINC01706242, and ZINC20552561) violated the Muegge rule, so these three were omitted for further analysis. The bioavailability of ligands was predicted using the swissADME tool [28]. Bioavailability is the concentration of a drug that can be absorbed by the body with respect to the systemic circulation and total dose. The bioavailability prediction score of a compound ranges from 0 to 1. Low bioavailability of a compound is indicated by the value of less than 0.5, whereas if the value is above 0.5, the compound is predicted to have high bioavailability. Based on the results of the predicted bioavailability, all analogues showed higher bioavailability (bioavailability score: 0.55). Therefore, all these 31 candidate analogues were selected for further study.

Predicted ADMET properties of these 31 analogues are shown in Table 2 . Positive HIA of the selected analogues indicate that these molecules are highly absorbed in the human intestine. The BBB is also positive, signifying that these compounds could enter the blood-brain barrier. CYPs play a vital role in the metabolism of drugs by catalyzing their biotransformation which is responsible for their degradation and excretion [52]; it is necessary to identify CYP-inhibiting properties during development of a new drug. Most of the analogues resulting from drug-likeness have low CYP450 inhibitory properties, except HCQ2, HCQ3, and CQ14 to CQ18 and CQ20 to CQ27. All the analogues except HCQ2, HCQ3, HCQ4, CQ15, CQ24 and CQ27 are non-inhibitors of P-glycoprotein, so these 6 analogues were not considered for further analysis in the study. When the genotoxic criteria of all the drugs were evaluated, they were found to be non-carcinogenic. Therefore, it was determined that only 13 analogues (1 HCQ and 12 CQ) were suitable candidates for the development of potential alternative compounds for HCQ and CQ against COVID-19 viral targets.

Table 2

Prediction of ADMET profile for HCQ and CQ analogues. The green colour indicates Blood-Brain Barrier (BBB) and Human Intestinal Absorption (HIA): Positive, P-glycoprotein Inhibitor (PGI): Non-inhibitor, CYP Inhibitory Promiscuity (CYP): Low, Carcinogens (C): Non-carcinogens and Aqueous solubility (AS): 0 to −4). The red colour indicates inhibitors of P-glycoprotein and Cytochrome P450 and high aqueous solubility.

The molecular docking studies showed that the spike glycoprotein interacts with HCQ1 through five amino acid residues, namely Pro426, Pro463, Phe464, Ser514, and Tyr396 with a binding energy of −6.83 kcal/mol (Table 3 , Fig. 1 ). In our study we found that the binding trend of the HCQ1 was higher than that of the parent drug; HCQ thus, this analogue displayed better binding affinity than that of the control drug. This analogue was found to interact with the binding pocket (Pro426, Pro463, Phe464, Ser514, Tyr396) as compared to that of the parent drug (Pro426, Phe464, Pro463, Glu516) suggesting similar interaction pattern. A pi-alkyl interaction occurs between HCQ1 and amino acid residues Pro426, Pro463 and Tyr396 of spike proteins. Conformational stability of the analogue HCQ1 is also associated with Pi-Pi T-shaped interactions with the amino acid residue, Phe464 and Pi-Donor Hydrogen Bond with amino acid residue, Ser514.

Table 3

Binding energy and the different types of interaction between the spike protein (PDB ID: 6LZG) of SAR-CoV-2 with CQ & HCQ and their analogues selected in the study.

Compounds	Binding Energy (kcal/mol)	Pi-Donor Hydrogen Bond	Pi-sigma	Pi-Alkyl Interaction	Pi-Pi T shaped	Pi-Pi Stacked	Pi-Anion	Van der Waals
HCQ	−3.98			Pro426, Phe464, Pro463, Glu516
HCQ1	−6.83	Ser514		Pro426, Pro463, Tyr396	Phe464
CQ	−4.38			Pro426, Asp428, Phe429, Pro463, Phe464
CQ1	−7.29		Thr430	Pro426, Asp428,Pro463, Tyr396	Phe46, Phe515			Glu516
CQ2	−6.67			Pro426, Asp428, Tyr396		Phe464	Glu516
CQ4	−6.07	Ser514		Pro426, Asp428, Phe464,Tyr396
CQ3	−6.05		Trp436	Phe338, Phe342,Phe374, Leu368, Val367	Phe464
CQ5	−6.03		Pro426	Asp428, Pro463	Phe464
CQ9	−5.89			Tyr396,Pro426, Asp428, Pro463			Glu516
CQ8	−5.26			Pro426, Asp428, Phe464,			Glu516
CQ6	−5.15	Phe338		Lue335, Cys336,Val362, Asp364
CQ12	−4.95			Asn343, Phe342, Leu368, Leu441		Phe374
CQ10	−4.84			Tyr396, Pro426, Pro463	Phe464
CQ20	−4.76		Ala372	Tyr369, Phe374
CQ7	−4.49		Trp436	Val367, Leu368, Phe374

Fig. 1

View of the top-ranking analogues (A1-B1) HCQ1 and (A2-B2) CQ1 docked in the ligand-binding site of the SAR-CoV-2 spike protein receptor. (A) Ligands with spike protein (hydrophobicity surface) at the active binding site. (B) 3D view of ligands with surrounding amino acids of the spike protein of SAR-CoV-2.

Similarly, all the 12 analogues of CQ show a better binding affinity than that of the parent drug, and their binding energies are in the range between −4.49 kcal/mol and −7.29 kcal/mol. Among the CQ analogues, CQ1 (−7.29 kcal/mol), CQ2 (−6.67 kcal/mol) and CQ4 (−6.07 kcal/mol) were the ranking inhibitors of SARS-CoV-2 spike glycoprotein. (Table 3, Supplementary Figure 1). From the results so obtained, it is clear that HCQ analogue such as “1-[1-(6-Chloroquinolin-4-yl) piperidin-4-yl] piperidin-3-ol” and CQ analogue- (1R,2R)-2-N-(7-Chloroquinolin-4-yl)cyclohexane-1,2-diamine shows higher binding affinity with the spike protein of the SAR-CoV-2. The interaction of these ligands with the receptor protein has been depicted in Fig. 1 and Table 3. CQ1 has the highest binding affinity and interacts with Pro426, Asp428, Pro463, Phe464, Phe515, Thr430, Glu516, Tyr396 amino acid residues of the spike glycoprotein as shown in Fig. 1. The stability of CQ1 at the interaction site of the spike protein is attributed to the presence of various Pi-interactions, such as Pi-alkyl interaction with amino acid residues, Pro426, Asp428, Pro463, Tyr396, Pi-Pi interaction with amino acid residues, Phe464, Phe515, and Pi-Sigma interaction with amino acid residue Thr430. MD analysis was carried out for the stable complexes of CQ and HCQ and their analogues for a period of 200 ns simulations run to understand the dynamic behaviour of the spike protein and to compare the difference between the parent and the analogue drugs. RMSD was plotted to understand the stability of the protein-drug complexes. The average RMSD values were observed to be 0.36 nm, 0.45 nm, 0.36 nm and 0.38 nm for Spike-CQ, Spike-HCQ, Spike-CQ1 and Spike-HCQ1, respectively. The RMSD plot (Fig. 2 ) suggests that the binding of the Spike-CQ, Spike-CQ1, Spike-HCQ remain stable throughout the entire MD simulation period. A sudden fluctuation is seen in Spike-HCQ1 complex at around 80 ns, but afterwards it is stable again till the end of the simulations. The CQ1-spike complex showed higher stability than that of the HCQ1 during the entire 200 ns simulation run.

Fig. 2

The root mean square deviation (RMSD) graph for the parent and analogue drug – spike protein complexes. Colour scheme: CQ, black; HCQ, red; CQ1, green; HCQ1, blue.

To gain an insight into the structural flexibility, root mean square fluctuation (RMSF) was examined. RMSF of the backbone atoms of the residues was plotted against simulated trajectory (Fig. 3 ). The average RMSF for all the complexes was evaluated. The HCQ1-Spike complex was observed to be 0.14 nm and the fluctuation was observed to be higher, with an average of 0.18 nm in the CQ1-Spike complex. The plot indicates global changes rather than local changes. Both RMSD and RMSF emphasize the stability aspects of the spike protein which seems to be impacted by the interaction of these selected drugs.

Fig. 3

The root mean square fluctuation (RMSF) graph for the drug-spike protein complexes. Colour scheme: CQ, black; HCQ, red; CQ1, green; HCQ1, blue.

The compactness of these complexes was also explored by analyzing their radius of gyration (Rg) values (Fig. 4 ). The average Rg values for CQ-Spike, HCQ-Spike, CQ1-Spike and HCQ1-Spike were observed to be 1.70 nm, 1.72 nm, 1.75 nm and 1.74 nm, respectively. The binding of the analogues showed higher Rg values than those of the parent drug. Lower Rg values were reported in case of CQ-Spike and HCQ-Spike complexes. A higher Rg value for selected analogues CQ1 and HCQ1was obtained in comparison to the parent CQ & HCQ-spike protein complexes indicating that the spike protein structure is considerably affected leading to the inhibition of cascading events. The solvent-accessible surface area (SASA) is characterized by the region of a protein accessible to solvent molecules. The average SASA values for CQ-Spike, HCQ-Spike, CQ1-Spike and HCQ1-Spike were 94.57 nm2, 95.51 nm2, 100.68 nm2, and 97.58 nm2, respectively during the 200 ns MD simulations. Results showed a comparable range in SASA values as observed with different spike complexes, particularly in Spike-CQ1 and Spike-HCQ1 complexes, which showed higher SASA values, whereas parent drugs CQ and HCQ with spike show lower SASA value (see Fig. 5).

Fig. 4

The radius of gyration (Rg) plot for the drug-spike protein complexes. Colour scheme: CQ, black; HCQ1, red; CQ1, green; HCQ1, blue.

Fig. 5

The Solvent Accessible Surface Area (SASA) plot for the drug-spike protein complexes as a function of time. Colour scheme: CQ, black; HCQ, red; CQ1, green; HCQ1, blue.

We also studied the collective motion of the analogues (HCQ1 and CQ1) and parent drugs (HCQ and CQ) with spike protein from their MD trajectories using principal component analysis (PCA). PCA was examined for understanding the structural and conformational changes in spike protein due to the binding of analogues and parent drugs. In this method, the dynamics of CQ-Spike, CQ1-Spike, HCQ-Spike and HCQ1-Spike were ascertained utilizing gmx covar module with reference to the backbone. The atomic fluctuations during eigenvector calculation were also reported in this study. The 2D projection of eigenvectors indicated diverse projections of spike protein upon binding with different molecules of ligands (Fig. 6 ). This analysis indicated the random fluctuations in the atoms of spike protein upon binding of different ligands and also suggests that the distinct positions in the spike protein undergo conformational changes.

Fig. 6

2D projection of trajectories on eigenvectors showed different projections of spike protein in case of (a) Spike-CQ (b) Spike-HCQ (c) Spike-CQ1 and (d) Spike-HCQ1. Colour scheme: CQ, black; HCQ, red; CQ1, green; HCQ1, blue.

To visualize the energy minima landscape of CQ, CQ1, HCQ and HCQ1 bound spike protein, we studied the free energy landscape (FEL) against first two principal components, Rg and RMSD which revealed ΔG value from 0 to 10 kJ/mol (Fig. 7 ). The shape and size of the minimal energy area (shown in blue) indicate the stability of the protein and protein-ligand complexes. Compact, unscattered and more centralized blue areas indicate the stability of the protein complex. Fig. 7 suggests that the analogue CQ1 bound spike complex is more stable than that of its parent drug, CQ-spike complex, and also than that HCQ and HCQ1 bound spike protein. Thus these analogues have the potential to induce spike protein to enter the local energy minimal state.

Fig. 7

Free Energy Landscape (FEL in kcal/mol) of the first two principal components for (a) Spike-CQ (b) Spike-HCQ (c) Spike-CQ1 and (d) Spike-HCQ1.

The basicity and acidity values for the listed top-ranking docked derivatives of HCQ and CQ including the parent compounds are falling in the range of 1.22–1.30 and 3.05 to 3.25, respectively (Table 4 ). This reflects that the analogues also exhibit similar biochemical response on viral inhibition by endosomal acidification to stop viral replication [45]. Thus, our study suggests that these selected analogues have a higher potential to work effectively as a novel antiviral analogue drug as compared to that of CQ and HCQ.

Table 4

Calculated acidity and basicity of CQ & HCQ and their analogues.

Compounds	Basicity	Acidity
HCQ	1.26	3.25
HCQ1	1.23	3.20
CQ	1.27	3.16
CQ1	1.24	3.31
CQ2	1.24	3.24
CQ4	1.26	3.23
CQ3	1.23	3.24
CQ5	1.30	3.04

Conclusions

Developing a robust framework for finding a desirable compound that binds and inhibits the attachment/internalization of SARS-CoV-2 might lead to potential therapeutics. This will overcome time constraints on antiviral drug discovery processes. Computational predictive approaches have been effectively used to develop a faster a prospective drugs or inhibitor candidates. In this study, we demonstrated an integrated computational screening framework to identify drug molecules that selectively target the spike protein of SARS-CoV-2 as a potent and safer alternative for CQ and HCQ. Virtual screening of similar compound libraries of CQ and HCQ was performed and the screened compounds along with parent drugs have been evaluated for drug likeness, ADMET properties, molecular docking and molecular dynamics simulations. Our docking study revealed that the identified analogue CQ1 has better binding affinity than that of HCQ1 and also that of parent drugs. Subsequently, MD simulations and PCA analysis of these complexes reform the stability of its mode of binding and thereby, inhibiting the activity of spike protein. Our findings suggest that the systematic selection of SARS-CoV-2 inhibitors was carried out by narrowing down on compounds and identifying from among them a suitable drug for COVID-19 treatment.

Author's contribution

RP, MS, and AG designed the study. PS, SSC, SP, MS, SG, and AG carried out the in silico work. MS, AG and SG wrote the manuscript, carried out data analysis for tabulation and prepared the figures. All the authors read and approved the final manuscript.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.