An in-silico evaluation of COVID-19 main protease with clinically approved drugs.

A novel strain of coronavirus, namely, SARS-CoV-2 identified in Wuhan city of China in December 2019, continues to spread at a rapid rate worldwide. There are no specific therapies available and investigations regarding the treatment of this disease are still lacking. In order to identify a novel potent inhibitor, we performed blind docking studies on the main virus protease Mpro with eight approved drugs belonging to four pharmacological classes such as: anti-malarial, anti-bacterial, anti-infective and anti-histamine. Among the eight studied compounds, Lymecycline and Mizolastine appear as potential inhibitors of this protease. When docked against Mpro crystal structure, these two compounds revealed a minimum binding energy of -8.87 and -8.71 kcal/mol with 168 and 256 binding modes detected in the binding substrate pocket, respectively. Further, to study the interaction mechanism and conformational dynamics of protein-ligand complexes, Molecular dynamic simulation and MM/PBSA binding free calculations were performed. Our results showed that both Lymecycline and Mizolastine bind in the active site. And exhibited good binding affinities towards target protein. Moreover, the ADMET analysis also indicated drug-likeness properties. Thus it is suggested that the identified compounds can inhibit Chymotrypsin-like protease (3CLpro) of SARS-CoV-2.

Introduction

The outbreak of the novel coronavirus disease, COVID-19, caused by coronavirus 2019–nCoV officially designated as severe acute respiratory syndrome related coronavirus 2 (SARS-CoV-2), represents a pandemic threat to global public health [1,2]. On January 30, the World Health Organization (WHO) announced a Public Health Emergency of International Concern (PHEIC) for the 2019-nCoV outbreak [3].

Coronaviruses are relatively large viruses containing a single stranded positive-sense RNA genome encapsulated within a membrane envelope. The viral membrane is studded with glycoprotein spikes that give coronaviruses their crown-like appearance. There are four classes of coronaviruses, namely alpha, beta, gamma and delta. The beta-coronavirus class include severe acute respiratory syndrome virus (SARS-CoV), the middle east respiratory syndrome virus (MERS-CoV), and the newly identified severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [4]. Although SARS-CoV-2 is classified into beta-coronaviruses group, it is diverse from MERS-CoV and SARS-CoV. It has been reported that SARS-CoV-2 genes share less than 80% identity (nucleotides) with SARS-CoV and it is more transmissible than other SARS-CoV viruses [[5], [6], [7]].

The SARS-CoV-2 genome consists of approximately 30,000 nucleotides. It encodes several structural proteins such as the glycosylated spike protein (S), envelope protein (E), membrane protein (M), and nucleocapsid protein (N). Additionally, the viral genome also encodes numerous nonstructural proteins, including RNA-dependant RNA polymerase (RdRp), coronavirus main protease (Mpro), and papain-like protease (PLpro). Upon entrance in the host cell, the viral genome is released and subsequently translated into viral poly-proteins using host cell translation machinery [1,6]. The poly-proteins are then cleaved into effector proteins by viral proteases PLpro and Mpro [8,9]. The Mpro, also known as 3-chymotrypsin-like protease (3 CL), plays a critical role in the virus replication process. It cleaves pp1a and pp1b polyproteins subsequently releasing functional proteins including RNA polymerase, endoribonuclease and exoribonuclease. Therefore, Mpro is a potential target for anti-coronaviruses screening. Indeed, inhibiting Mpro activity could stop the spread of infection [10,11].

The released crystal structure of Mpro (6lu7) was obtained by a co-crystallization with a peptide-like inhibitor called N3 (PRD_002214). The enzyme has a molecular weight of 33.79 kDa as determined by Mass spectroscopy. Mpro forms a dimer, where each monomer comprises three domains: Domain I (residues 8–101) and II (102–184) consists of an antiparallel beta barrel, and the alpha helical domain III (residues 201–301) is required for enzymatic activity. It shares a similar structure with cysteine protease, however the active site lacks the third catalytic residue. Thus, catalytic harbors Histidine 41 (H41) and Cysteine 145 (C145) [12,13] located in between domain I and II, while amino acids, T24, L27, H41, F140, C145, H163, M165, P168 and H172 form a hydrophobic surrounding in the pocket [14].

The concept of drug repurposing is widely used nowadays to identify potential drugs for COVID-19 disease. The practice of reusing drugs for more than one purpose can significantly reduce the cost, time and risks of the drug development process [15].

To date many studies have targeted SARS-CoV-2 Mpro, out of which, the study of Khan and collaborators, based on computational drug design methods, yielded five promising hits corresponding to three antiviral drugs (Remdesivir, Saquinavir and Darunavir) and two natural compounds (Flavone and coumarine) [16]. Similarly, Aanouz and collaborators proposed three compounds (Corcin, Digitoxigenin, and β-Eudesmol) as potential inhibitors of Mpro from Moroccan medicinal plants [17]. Interestingly, Talampicillin and Lurasidone are reported to possess the highest affinity towards Mpro enzyme and might be repurposed against COVID-19 [11].

Although, several molecular docking studies have been established to find a potential inhibitor of Mpro activity based on antiviral compounds commonly used to treat immunodeficiency virus (HIV), phytochemical, antimalarial agents, spices or marine products [[18], [19], [20], [21], [22]]. Unfortunately, no specific therapies for COVID-19 disease are available and investigations regarding the treatment of COVID-19 disease are still lacking. Therefore, there is a dire need to identify approved drugs which may inhibit SARS-CoV-2 virus proteins and to find at the same time their optimal association and concentration to treat COVID-19 patients.

In the following study, we investigated the Mpro inhibitory potential of eight clinically approved drugs belonging to four pharmacological classes: anti-viral, anti-bacterial, anti-infective and anti-histamine, using a molecular docking and a molecular dynamics approach.

Methods

Target and ligand retrieval form databases

This research is a descriptive-analytical study. In this study, the interaction of several approved compounds was investigated. A total of eight compounds were tested against COVID-19 main protease (Mpro). N3 compound was used as a docking target for comparison.

In order to obtain the structure information of selected compounds, a Drugbank database (https://www.drugbank.ca/) was used [23] (Table 1 ). The Three dimensional structure of Mpro (Fig. 1 ) was retrieved from the Protein Data Bank (PDB) (https://www.rcsb.org/) [24]. It corresponds to a complex between the enzyme and its inhibitor N3 [13]. The 6lu7 structure preparation consists of several steps such as deleting all water molecules, N3 inhibitor and adding hydrogens. The new file was saved for docking analysis.

Table 1

List of approved compounds used in this study.

Compound	Accession numbera	Drug class	MM g/mol	Description
Chloroquine	DB00608	Aminoquinolone derivative	319.87	Antimalarial
Quinine	DB00468	Alkaloid	324.41	Antimalarial
Nitazoxanide	DB00507	Thiazolide	307.28	Anti-infective
Doxycycline	DB00254	Oxytertacycline derivative ATBb	444.43	Antibacterial
Lymecycline	DB00256	Tetracycline with a 7-chloro substitution	602.63	Antibacterial
Cetirizine	DB00341	Histamine H1 antagonist	388.88	Antihistamine
Mizolastine	DB12523	Histamine H1 antagonist	432.50	Antihistamine
Indinavir	DB00224	Specific HIV protease inhibitor	613.78	Antiviral

Drug bank accession number.

Antibiotic, MM: Molecular Mass.

Fig. 1

Pymol Structure representation of Mpro. Surface (A) and Cartoon (B) representations of one protomer of the dimeric CoV Mpro. (C) Surface and stick representations of the conserved binding pocket of Mpro. Red color: Mpro domain I, Bleu color: Mpro domain II and pink color for Mpro domain III. Gray color represents coils. C(Cys), D (Asp), E (Glu), F(Phe), H(His), L (Leu), M(Met), Q (Gln). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Blind docking

In order to investigate the molecular interaction between several approved compounds and COVID-19 main protease (Mpro), blind docking was performed using a SwissDock server (http://www.swissdock.ch/) under the accurate mode with no flexibility of the side chain of any amino acid of the target protein. In addition, a binding pocket was not defined so as not to bias the docking towards the active site. The protein and ligand were specified by uploading PDB and Mol2 files, respectively.

SwissDock generates all possible binding modes for each ligand and the most favorable binding modes at a given pocket were clustered. All ligand clusters were saved in an output file called “prediction file”. The prediction file provided; Cluster rank, Element, Fulfitness and estimated binding free energy (ΔG). A cluster corresponds to a predicted binding pocket on the target protein and the cluster rank represents the different conformations of a ligand in a certain cluster [25]. Only the lowest energy model of cluster zero was considered to be the most favorable interaction.

Clusters visualization and 2D diagrams generation

After docking Chimera and Pymol softwares were used to visualize the receptor ligand interactions for the lowest energy model of the clusters obtained from the previous step. Each ligand cluster was inspected for amino acids interacting with the ligand, hydrogen bonds (H-bonds), the specific atoms involved. All the interacting amino acids with the target were noted for each cluster [26,27].

2D diagrams were generated by PoseView integrated into Proteins Plus server (https://proteinsplus.zbh.uni-hamburg.de/#dogsite) [28], which automatically creates two-dimensional diagrams of protein ligand complexes according to chemical drawing conventions. The generation of structure diagrams and their layout modifications are based on the library 2Ddraw. Interactions between the molecules are estimated by a builtin interaction model that is based on atom types and simple geometric criteria [29,30].

Binding pocket prediction

Although the binding site is well characterized for N3 within many CoV Mpro crystals. We have applied the DoGSiteScorer online tool (https://proteinsplus.zbh.uni-hamburg.de/#dogsite) to predict and describe binding pockets within native Mpro and the complexes Mpro/inhibitors obtained after docking analysis.

DoGSiteScorer is a grid-based method which uses a Difference of Gaussian filter to detect potential binding pockets solely based on the 3D structure of the protein and splits them into subpockets. Global properties, describing the size, shape and chemical features of the predicted (sub) pockets are calculated. Per default, a simple druggability score is provided for each (sub) pocket, based on a linear combination of the three descriptors describing volume, surface, hydrophobicity and enclosure. Furthermore, a subset of meaningful descriptors is incorporated in a support vector machine (libsvm) to predict the (sub) pocket druggability score (values are between zero and one). The binding pockets are ranked according to their size, surface area and druggability score [31].

Drug likeliness, pharmacokinetic and oral toxicity evaluations

To analyze drug likeliness and the pharmacokinetics parameters, all the molecules were subjected to ADMET (Absorption, distribution, metabolism, excretion and toxicity) predictions with SwissADME (http://www.swissadme.ch/) and PROTOX web servers (http://tox.charite.de/protox_II/) [32,33]. PROTOX predict the median oral lethal dose, LD50 (Lethal dose, 50%) based on the analysis of the two-dimensional similarity to compounds with known LD50 values and the identification of fragments over-represented in toxic compounds [33].

Molecular dynamics simulation

To gain atomic level insight into the binding interaction of Lymecycline and Mizolastine with SARS-Cov 2 Mpro, all atom-molecular dynamics simulations were carried out using PMEMD. cuda module implemented in AMBER 18 with explicit water model under periodic boundary conditions [34]. The parameters for both systems (Mpro-Lymecycline and Mpro-Mizolastine) were prepared using the antechamber and tleap modules. Amber ff14SB [35] force field was assigned to protein whilst Amber Generalised Force Field (GAFF) [36] was selected for the two compounds. Each system was solvated using the TIP3P water model in a periodic box, with a minimum spacing of 10 Å from the solute/protein atoms. After that neutralization of the system with counter ions (Na+) replacing solvent molecules at the position of electrostatic favorable potential. The Particle Mesh Ewald Method (PMEM) and the SHAKE algorithm were used to calculate the long-range electrostatic interactions and to constrain the hydrogen bonds respectively [37,38]. To remove bad contacts and possible steric clashes both complexes were subjected to energy minimization. Briefly systems were relaxed by adjusting hydrogen position, whereas for the Mizolastine complex the first 2500 cycles were carried out using steepest descent algorithm and the next 2500 steps involved the conjugate gradient algorithm. However, for Lymecycline complex the number of cycles were increased to achieve minimized system. The minimization steps were repeated multiple times with gradual decrease in force applied to restrain heavy atoms. It was followed by 5000 (50,000 for lymecycline complex) additional minimization steps in absence of any restraint. Afterwards each system was equilibrated under NVT (gradual heating from 0 K to 300 K) and NPT (pressure = 1 atm) ensembles. The final production run of 120 ns was performed at constant temperature (300 K) and pressure (1 atm) with integration time steps of 2 fs. The obtained results were analyzed via CPPTRAJ modules. Chimera and VMD were used for visualization and the graphs were plotted using Xmgrace tool [27,39].

Energy calculations

To estimate the binding free energies of the selected inhibitors against SARS-CoV-2 Mpro, the widely known MM-PBSA was performed. This tool computes components of binding free energy utilizing the molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA), implemented in AMBER 18. A total of 2500 frames were extracted from the most stable trajectories and of them 500 frames were used to compute various energy components involved in ligand binding.

In the MMPBSA protocol, the following equation is used to calculate binding free energy (ΔGbind) of an inhibitor.ΔGbind = ΔH−TΔS ≈ΔEMM + ΔGsolv − TΔS ΔEMM = ΔEinternal + ΔEelec + ΔEvdW ΔGsolv = ΔGpol + ΔGnpWhere ΔEMM represents the alterations in molecular mechanics energy at gas phase, ΔGsolv and -TΔS indicate desolvation free energy and the conformational entropic contributions upon binding with ligand, respectively. Further, ΔEMM is sum of ΔEinternal (bond, dihedral and angle energies), ΔEelec (electrostatic), and ΔEvdW (Van der Waals) energies while ΔGpol and ΔGnp represents polar and nonpolar solvation energy [40].

Results and discussion

Molecular docking

In order to investigate the possible mechanism by which selected drugs act, an in silico theoretical molecular docking approach was used. Fig. 2 illustrates docking poses of the studied compounds.

Fig. 2

Docking poses of different drugs against Mpro visualized by Pymol. The protease Mpro is shown as gray background, inhibitors are in different colors. (A) Indinavir. (B) Chloroquine. (C) Lymecycline. (D) Mizolastine. (E) Quinine. (F) Cetirizine. (G) Nitazoxanide. (H) Doxycycline. H-bonds are represented by black dashed lines. Interacting residues are labeled: E (Glu), G (Gly), H (His), L (Leu), N (Asn), Q (Gln), T (Thr). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

During our study, we simulated the binding mode of N3 against 6lu7 crystal structure using SwissDock to ensure the effectiveness of docking results and to compare results produced by several drugs to those of N3. Indeed, this compound is a well characterized inhibitor of COVID-19 main protease.

Docking results revealed that N3, Indinavir and Chloroquine had the best energies of binding −10.83, −9.81 and −9.71 kcal/mol, respectively (Table 2 , column 5), which is consistent with three studies. The first one reported the entire complex N3/Mpro crystal structure saved in the PDB database under 6lu7 accession number [13]. The second reported that Indinavir exhibited a good docking score (−7.05) when docked against 5r7z Mpro structure using flexible docking with Glide and the last one revealed that Chloroquine and its derivatives can bind to Mpro [[18], [21]].

Table 2

Molecular docking analysis results for several drugs against 6lu7 crystal structure. These drugs were ranked according to their minimum binding energy. The lowest energy model of cluster rank zero was considered.

Compounds	Clusters	Total elements	Fullfitness (kcal/mol)	ΔG (kcal/mol)	H-bonds	Length (Å)	Residues
N3	38/38	256/256	−1172.91	−10.83	2	2.1	G143
						2.8	Q189
Indinavir	54/54	256/256	−1098.93	−9.81	2	2.7	L141
						2.6	G143
Chloroquine	44/44	254/254	−1223.94	−9.71	1	2.6	N142
Lymecycline	23/33	168/256	−1332.56	−8.87	3	2.2	E166
						2.3	E166
						2.9	Q189
Mizolastine	54/54	256/256	−1300.12	−8.71	3	2.3	G143
						2.5	T24
						2.1
Quinine	39/40	242/250	−1135.13	−8.09	2	2.6	E166
Cetirizine	38/42	224/256	−1112.62	−7.99	1	2.1	N142
Nitazoxanide	58/64	224/256	−1215.20	−7.77	1	2.1	E166
Doxycycline	22/32	176/256	−1276.83	−7.52	2	2/2.4	E166

Column 2 represents clusters within Mpro binding pocket/total clusters.

Column 3 represents the number of conformations within the Mpro binding pocket/Total elements. E(Glu), G(Gly), L(Leu), N(Asn), Q(Gln), T(Thr).

According to a full fitness score, Lymecycline and Mizolastine had more favorable binding mode, indicated by more negative fullfitness scores −1332.56 and −1300.12 kcal/mol, respectively (Table 2, column 4). Azelastine can undergo an optimal binding with Mpro [41]. When the clusters were analyzed it was found that N3, Indinavir, Chloroquine and Mizolastine showed that all clusters were able to fit into Mpro binding pocket. However, Lymecycline revealed a binding energy of −8.87 kcal/mol and it was able to occupy 23 clusters constituting a total of 168 possible conformations within the substrate binding cavity, out of a total of 256 elements (Table 2, column 3). It is striking that all or more than half of the total predicted elements are docked in the substrate binding pocket.

To investigate the possible reasons for differences in the binding energies, we examined the docked complexes with Pymol software and PoseView integrated in Protein Plus server. Table 2 showed the number and length of H-bonds formed between the target protease and the different compounds. Chloroquine, Nitazoxanide and Cetirizine established only one H-bond with N142, E166 and N142 residues, respectively. Otherwise, Quinine, Doxycycline and Indinavir were found to form two H-bonds with E166, E166 and (L141, Gly 143) residues, respectively. Interestingly, Mizolastine and Lymecycline were found to form three H-bonds with Mpro. Indeed, Lymecycline established three specific H-bounds (two H-bonds with E166 residue and one H-bond with Q189 residue) (Fig. 3 A, and C). However, Mizolastine formed three H-bonds (two H-bonds with T24 residue and one with G143) and hydrophobic interactions with N142 and C145 residues (Fig. 3B and D).

Fig. 3

Interactions established after docking of Lymecycline and Mizolastine with 6lu7 Mpro protease. Three dimensional representation of intermolecular interaction in (A) Complex Lymecycline/6lu7 Mpro, (B) Complex Mizolastine/6lu7 Mpro (Cyan dashed lines depict hydrogen bonds). 2D plot showing both hydrogen and hydrophobic interactions in (B) Complex Lymecycline/6lu7 Mpro, (C) Complex Mizolastine/6lu7 Mpro. Directed bonds between protein and ligand are drawn as black dashed lines and the interacting protein residues and the ligand are visualized as structure diagrams. Hydrophobic contacts are represented more indirectly by means of spline sections highlighting the hydrophobic parts of the ligand and the label of the contacting amino acid. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

N3 compound forms multiple hydrogen bonds with the main chain of residues in the substrate-binding pocket [31]. However, only two hydrogen bonds were detected after SwissDock analysis because the complex of Mpro with its inhibitor N3 was obtained in theoretical (in silico) not in experimental conditions.

Although Lymecycline and Doxycycline belong to the same family of tetracyclines, Lymecycline bind more effectively to Mpro with a minimum energy of −8.87 kcal/mol compared to Doxycycline with −7.52 kcal/mol binding energy, suggesting, the role of NH (CH2)4 CH COOH NH2 chemical substituting group in increasing the binding affinity of Lymecycline towards Mpro enzyme.

Binding pockets prediction

To elucidate and describe binding pockets within Mpro, Dogsitescorer server was used to analyze unliganded structure (6lu7) and N3/6lu7, Lymecycline/6lu7, Mizolastine/6lu7 complexes. Here, we have reported only the first three pockets in Table 3 . Results revealed that N3, Lymecycline and Mizolastione occupied the same pocket (P0) within 6lu7 structure (Fig. 4 ) with a high druggability score of 0.8 and a volume of 1191.74, 1061.18 and 1266.18 Å3, respectively.

Table 3

Binding pockets prediction for chain A of Mpro structure (6lu7) in complex with inhibitors (N3, Lymecycline and Mizolastine).

Structure	Pocket Number	Volume (Å3)	Surface (Å2)	Drug score	Hydrophobicity ratio	Enclosure
6lu7	0	702.27	842.81	0.77	0.35	0.13
	1	374.59	757.16	0.74	0.48	0.11
	2	330.18	518.79	0.56	0.40	0.24
6lu7+N3	0	1191.74	1136.13	0.8	0.34	0.1
	1	257.79	538.69	0.56	0.43	0.08
	2	253.7	544.53	0.51	0.49	0.24
6lu7+Lymecycline	0	1061.18	1032.51	0.8	0.35	0.08
	1	292.67	475.29	0.51	0.39	0.25
	2	277.7	601.59	0.65	0.47	0.1
6lu7/Mizolastine	0	1266.18	1294.39	0.8	0.32	0.09
	1	272.13	514.16	0.64	0.34	0.1
	2	272.06	488.48	0.5	0.49	0.24

Fig. 4

Pymol solid surface representations of binding mode of N3 compound (A), Lymecycline (B) and Mizolastine (C). Binding pocket prediction of N3 (D), Lymecycline (E) and Mizolastine (F) within 6lu7 pocket. Amino acid residues are labeled: N (Gln) and H (His).

Based on Docking and Dogsitescorer studies, it is evident that Lymecycline and Mizolastine showed favorable binding with the new Mpro, and the data comparable with those of N3.

Pharmacokinetics properties

Drug likeness and pharmacokinetics proprieties analysis showed that Mizolastine and Chloroquine had four interesting features; a high lipophilicity, bioavailability score of 0.55, and higher gastrointestinal (GI) absorption. Otherwise, Lymecycline is soluble compared to Chloroquine, Indinavir and N3. The absorption is low and presents a bioavailability score of ≈0.1 with N3 compound. Interestingly, Lymecycline and Chloroquine are not substrates of P- glycoprotein compared to other drugs (Table 4 ). Indeed, P-glycoprotein is a multidrug transporter which is apically expressed in the gastrointestinal tract, liver, kidney and brain endothelium. Consequently, P-glycoprotein plays an important role in the oral bioavailability, CNS (Central Nervous System) distribution biliary and renal elimination of drugs which are substrates of this transporter. It has been reported that Indinavir and Saquinavir displayed a reduced antiviral activity. Both drugs partially recover the ability to inhibit the replication of HIV-1 in the presence of Verapamil (P-glycoprotein Inhibitor). Thus, P-gp expression may affect the activity of antiviral drugs [42,43].

Table 4

Drug likeliness, pharmacokinetic and LD50 proprieties of the eight selected drugs

P-gp: P-glycoprotein, GI: Gastrointestinal, BBB: Blood Brain Barrier, LD: Lethal dose, 50%.

Compounds	Water solubility	Lipophilicity Consensus Log Po/w	GI absorption	BBB permeability	P-gp substrate	Bioavailability score	LD50 mg/kg
N3	Poorly soluble	2.69	Low	No	Yes	0.17	4000
Indinavir	Poorly soluble	2.78	High	No	Yes	0.55	5000
Chloroquine	Poorly soluble	4.15	High	Yes	No	0.55	311
Lymecycline	Soluble	−1.86	Low	No	No	0.11	3000
Mizolastine	Poorly soluble	3.28	High	Yes	Yes	0.55	450
Quinine	Moderately soluble	2.81	High	Yes	No	0.55	263
Cetirizine	Moderately soluble	2.56	High	Yes	Yes	0.55	365
Nitazoxanide	Soluble	1.16	Low	No	No	0.55	1350
Doxycycline	Soluble	−0.34	Low	No	Yes	0.11	2240

Toxicity data analysis showed that Lymecycline (LD50 = 3000 mg/kg) and Mizolastine (LD50 = 450 mg/kg) are safe to used compared to Chloroquine (LD50 = 311 mg/kg). In addition, several studies reported that Tetracyclines had clinical utility against intracellular parasites such as malaria, Chlamydia and a-proteobacteria [44].

MD simulation analysis

Molecular dynamic is a state-of-the-art simulation method for studying the physical motion and trajectory of the atoms in the presence of other molecules along with the various interactions within a system. It helps to follow and understand the structural features and conformational dynamics in the system. Thus, to validate the stability of the system and to probe ligand induced perturbations, MD simulation was performed with two best compounds as a function of time. The MD trajectories were examined based on various parameters including Root Mean Square Deviation (RMSD), Root Mean Square Fluctuation (RMSF), Radius of Gyration (Rg), Inter-molecular hydrogen bond interaction and occupancy. Moreover, binding free energy calculations were also performed.

RMSD monitors the deviations in average distance between the atoms of target protein during simulation with respect to initial docking structure/reference frame. In short it is the deviation of the 3D structure over time. It provides insight into the system’s stability, equilibrium and convergence whereas, the smaller fluctuations and constant backbone atoms (C, Cα, N, and O) RMSD, is indicative of the stable system. As described in Fig. 5 A after an initial period of fluctuation both systems attained equilibrium during the last 50 ns of the simulation run. In general, the Mpro and Lymecycline system displayed slightly higher fluctuation, whilst in comparison the lowest deviations were observed for Mpro-Mizolastine complex. For Lymecycline complex during the initial frames continuous increase in RMSD value was observed in the range of <2 - 4 Å however, in the last 50 ns trajectories the system obtained stability with the deviation of <3 Å whereas no sharp fluctuations were observed during this time frame. In comparison, for Mizolastine complex, after gradual increase in fluctuation during the initial 45 ns time period, the system attained equilibrium state in the remainder MD trajectories except the frames in between 60 and 65ns where sharp fluctuation peaks yet in acceptable range (<3.8 Å) were observed. The average RMSD for Mpro-Lymecycline and Mpro-Mizolastine complex was maintained at 3.10 ± 0.43 and 3.66 ± 1.77 Å, which implies that both systems attained a more stable structure compared to initial structure. Additionally, there was not much deviation between average and observed RMSD of protein at the end of the 120 ns simulation and for both systems the RMSD during the whole run was <4 Å which is in acceptable range.

Fig. 5

Time evolution plots of Molecular Dynamics Simulation trajectories of Mpro-Lymecycline complex and Mpro-Mizolastine complex (A) Root Mean Square Deviation (RMSD), (B) Root Mean Square Fluctuation (RMSF) and (C) Radius of Gyration (Rg).

vRMSF analysis is highly suitable to calculate the time dependent fluctuations in each residue, in turn is an essential parameter to determine the protein’s flexibility. The RMSF of backbone atoms was calculated for Mpro with 306 residues with the two potential inhibitors. As illustrated in Fig. 5B, both complexes had almost similar fluctuations trend. For most of the protein residues the observed RMSF value was lower than 2.5 Å. Whereas, the residues interacting with the ligands in the active site were found stable and displayed little fluctuations over time indicating the stable nature of compounds with target protein. In contrast the main fluctuations corresponded to the region that were distant from the ligand active site and others were found around the flexible loop region. Furthermore, the RMSF fluctuations in the Mizolastine complex was observed to be lower than the Lymecycline complex, suggesting that it had relatively lower structural mobility than the Lymecycline. The average RMSF values were 1.46 ± 0.93 and 1.32 ± 0.57 for Lymecycline and Mizolastine respectively.

The Rg is a parameter to assess the compactness and overall dimension of the protein which in turn signifies folding and unfolding of the protein. The lower the gyration values the more folded the protein is and vice versa. Therefore, Rg was calculated to determine whether the tested drugs maintained the folding of the system. The gyration graph for backbone carbon atoms of the protein relative to time is presented in Fig. 5C. The Rg values for both the systems range from 22 to 22.5 Å with an average score of 22.31 ± 0.14 Å for Lymecycline and 22.43 ± 0.28 Å for Mizolastine. As evident from the plot, the Lymecycline complex attained more compacted form as the MD simulation progressed, indicating a well converged system.

Free energy calculations for protein-ligand binding

Free energy calculations are a collection of methods used to envisage ligand binding affinities by considering several atomic level interactions that could be responsible for the affinity of ligand toward targeted protein. Thus, to understand the biophysical basis of interaction between selected hits and main protease in further detail MM/PBSA calculations were carried out. For each complex non-polar solvation energy (Gnonpolar), polar solvation (GPB), electrostatic (Eelec), Van der waals (Evdw) and binding free energy (ΔGbinding) were calculated. A summary of individual components involved in the binding free energy of tested inhibitors with main protease is shown graphically in Fig. 6 with the data listed in Table 5 .

Fig. 6

The energy contributions (KJ/mol) of Lymecycline and Mizolastine in complex with SARS-CoV-2 main-protease (PDB ID 6lu7) in MM/PBSA assay.

Table 5

The binding energies (KJ/mol) for Lymecycline and Mizolastine in complex with SARS-CoV-2 main-protease (PDB ID 6lu7) using MM/PBSA from the snapshots extracted from MD simulation trajectories.

Energy Components	Complexes
	6lu7-Lymecycline	6lu7-Mizolastine
Evdw	−34.97 ± 3.85	−43.73 ± 3.24
Eelec	−93.63 ± 20.92	−30.16 ± 4.91
EPB	110.27 ± 18.99	50.34 ± 5.08
Enonpolar	−3.87 ± 0.25	−4.33 ± 0.20
ΔGgas	−128.60 ± 21.54	−73.89 ± 5.92
ΔGsolv	106.40 ± 18.87	46.02 ± 4.99
ΔTotal	−22.19 ± 5.23	−27.87 ± 3.91

It is evident from the data that the binding of Lymecycline to the target protein is favored mostly by intermolecular electrostatic interactions which is also evident from the numbers of hydrogen bond contacts formed between inhibitor and targeted enzyme. Van der Waals interactions and non-polar solvation energy also contributed in binding affinity. In contrast, for Mizolastine Van der Waals interactions were mainly responsible for complexation along with electrostatic and non-polar interactions. However, the complexation in each system is disfavored by polar solvation energy. The calculated ΔG binding energy for Lymecycline and Mizolastine with target protein was found to be −22.19 ± 5.23 and −27.87 ± 3.91 kcal/mol, respectively.

Inter-molecular interactions and hydrogen bond occupancy

Generally, Protein-ligand interaction is stabilized by different types of weak interactions with hydrogen bonding contacts being the most important interaction. Hydrogen bond analysis was performed for both the complexes, with the occupancies reported in Table 6 . In general H-bond patterns were observed with various active site residues during simulation. The inter-molecular interactions were plotted utilizing the structural ensembles from stable trajectories. In the Mpro-Lymecycline complex the active residues engaged in H-bond by Lymecycline were L141, G143, S144, catalytic dyad residue C145, E166. The highest occupancy was observed for E166 (36.1%), followed by S144 (34.7%). The molecular interaction diagram (Fig. 7 A) indicated that carbonyl oxygen on cyclohexene ring from Lymecycline contributed to the formation of hydrogen bonds with S144, and G143 amide hydrogen and also with the carbonyl oxygen of S144. Moreover, E166 displayed two H-bond contacts with the side chain of Lymecycline. Other H-bonds were observed between carbonyl oxygen of methylated cyclohexane ring from Lymecycline and amide hydrogen of C145 and with the L141 residue.

Table 6

Hydrogen bonds between Lymecycline and Mizolastine with SARS-CoV-2 Mpro that found with at least 0.5% occupancy throughout 120 ns simulation run.

Complexes	Hydrogen bond Formation Acceptor-DonorH/Donar	Distance (Å)	Occupancy (%)
Mpro-Lymecycline	GLU_166@O-Lymecycline@H15/Lymecycline@N14	2.80	36.1
	Lymecycline@O35-SER_144@HG/SER_144@OG	2.71	34.7
	Lymecycline@O35-GLY_143@H/GLY_143@N	2.87	24.8
	Lymecycline@O11-GLU_166@H/GLU_166@N	2.87	21.2
	Lymecycline@O35-SER_144@H/SER_144@N	2.89	19.1
	Lymecycline@O38-CYS_145@H/CYS_145@N	2.90	11.6
	GLU_166@OE1-Lymecycline@H25/lymecycline@N23	2.76	11.2
	Lymecycline@O38-GLY_143@H/GLY_143@N	2.83	11.0
	LEU_141@O-Lymecycline@H71/Lymecycline@O33	2.80	10.8
	GLU_166@OE1-Lymecycline@H26/Lymecycline@N23	2.75	8.85
	GLU_166@OE1-Lymecycline@H24/Lymecycline@N23	2.76	8.47
	Lymecycline@O28-ASN_142@H/ASN_142@N	2.85	6.82
	GLU_166@OE2- Lymecycline@H25/Lymecycline@N23	2.76	5.17
	GLU_166@OE1-Lymecycline@H3/Lymecycline@N2	2.79	5.16
	GLU_166@O-Lymecycline@H16/Lymecycline@N14	2.79	4.71
	PRO_168@O-Lymecycline@H79/Lymecycline@O48	2.74	4.57
	GLU_166@OE2-Lymecycline@H24/LIG_306@N23	2.76	4.32
	ASN_142@OD1-Lymecycline@H24/Lymecycline@N23	2.78	3.77
	GLU_166@OE2-Lymecycline@H26/Lymecycline@N23	2.77	3.55
	ASN_142@OD1-Lymecycline@H26/Lymecycline@N23	2.79	3.55
	ASN_142@OD1- Lymecycline@H25/Lymecycline@N23	2.79	3.29
	GLU_166@OE1-Lymecycline@H16/Lymecycline@N14	2.78	2.97
	GLU_166@OE2-Lymecycline@H3/Lymecycline@N2	2.81	2.40
	GLU_166@OE2-Lymecycline@H16/Lymecycline@N14	2.78	2.11
	Lymecycline@O29-ASN_142@H/ASN_142@N	2.85	1.85
	Lymecycline@O33-ASN_142@HD22/ASN_142@ND2	2.88	1.72
	ASN_142@OD1-Lymecycline@H71/Lymecycline@O33	2.71	1.46
	SER_144@OG-Lymecycline@H71/Lymecycline@O33	2.82	1.10
	GLU_166@OE2-Lymecycline@H15/Lymecycline@N14	2.79	1.09
	Lymecycline@O29-SER_46@HG/SER_46@OG	2.68	1.08
	PHE_140@O-Lymecycline@H25/Lymecycline@N23	2.83	1.0
	Lymecycline@O11-GLN_189@HE21/GLN_189@NE2	2.85	0.89
	THR_190@O-Lymecycline@H26/Lymecycline@N23	2.79	0.79
	GLU_47@OE1-Lymecycline@H25/Lymecycline@N23	2.75	0.69
	PRO_168@O-Lymecycline@H24/Lymecycline@N23	2.81	0.65
	GLN_189@OE1-Lymecycline@H16/Lymecycline@N14	2.82	0.55
	THR_190@O-Lymecycline@H25/Lymecycline@N23	2.79	0.52
Mpro-Mizolastine	Mizolastine@N15-GLY_143@H/GLY_143@N	2.92	18.8
	Mizolastine@O8-THR_24@HG1/THR_24@OG1	2.73	11.2
	Mizolastine@O8-SER_46@H/SER_46@N	2.86	10.0
	THR_24@O-Mizolastine@H38/Mizolastine@N9	2.88	3.54
	SER_46@OG-Mizolastine@H38/Mizolastine@N9	2.90	2.83
	Mizolastine@N15-ASN_142@HD22/ASN_142@ND2	2.92	1.17
	Mizolastine@O8-THR_45@HG1/THR_45@OG1	2.74	0.89
	Mizolastine@N4-THR_26@H/THR_26@N	2.95	0.83
	Mizolastine@N15-CYS_145@H/CYS_145@N	2.95	0.78
	Mizolastine@O8-THR_26@H/THR_26@N	2.86	0.57

Fig. 7

MD intermolecular interaction analysis. (A) Binding orientation of Lymecycline and mode of interaction with SARS-CoV-2 Mpro. (B) Binding orientation and mode interaction of Mizolastine with SARS-CoV-2 Mpro.

Mizolastine was found to make H-bond interaction with G143 (18.8%), N142 (1.17%) and C145 (0.78%) and T190 (0.30%) from the active site. Furthermore, T24, T26, T45 and S46 were also involved. To investigate the stability of interactions and to explore binding mode, interactions analysis was also carried out utilizing the same information used to calculate MMPBSA. As depicted in Fig. 7B three hydrogen bonds were observed. The Imidazole ring from benzimidazole scaffold contributes to the hydrogen bond with C145 and G143. Another H-bond was formed between the pyrimidinone ring and T26 (0.83%). It is important to note here that from the H-bond interactions analyzed in docking simulation only the interaction with G143 was maintained till the end of the simulation. For T24 the H-bond occupancy found was 11.2 and 3.45% however, the distance increased from 2.1 before MD simulation to 3.6 Å in the most stable structure during simulation.

Conclusion

This blind molecular docking study proposes potential available approved drugs: Lymecycline and Mizolastine as prospective inhibitors of SARS-CoV-2 Mpro. The estimated free energy of binding of these two drugs is in the range of −8.87 and −8.71 kcal/mol, proffering the spontaneous and energetically favored production of protein ligand complex. Molecular dynamic study concludes that both of the tested compounds, can serve as potential Mpro inhibitors based on their stable nature during a long term molecular dynamic simulation. Though the maximum number of hydrogen bonds were formed by Lymecycline during the simulation the Mizolastine was stabilized in the pocket mainly via non-bonded interactions as indicated from energy component analysis. In summary, the data from RMSD, RMSF, Rg and MMPBSA indicated that there was very insignificant/minor difference between the two in imparting stability to the system. Therefore, these are potential candidates for further in vitro testing and Anti SARS-CoV-2 therapeutics development.

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.