key: cord-0861242-mhibzr4d authors: Alamri, Mubarak A.; Tahir ul Qamar, Muhammad; Afzal, Obaid; Alabbas, Alhumaidi B.; Riadi, Yassine; Alqahtani, Safar M. title: Discovery of anti-MERS-CoV small covalent inhibitors through pharmacophore modeling, covalent docking and molecular dynamics simulation date: 2021-05-15 journal: J Mol Liq DOI: 10.1016/j.molliq.2021.115699 sha: 9292e94117b227cd52e0c6e4af396a37618f0b7f doc_id: 861242 cord_uid: mhibzr4d Middle east respiratory syndrome coronavirus (MERS-CoV) is a fatal pathogen that poses a serious health risk worldwide and especially in the middle east countries. Targeting the MERS-CoV 3-chymotrypsin-like cysteine protease (3CL(pro)) with small covalent inhibitors is a significant approach to inhibit replication of the virus. The present work includes generating a pharmacophore model based on the X-ray crystal structures of MERS-CoV 3CL(pro) in complex with two covalently bound inhibitors. In silico screening of covalent chemical database having 31,642 compounds led to the identification of 378 compounds that fulfils the pharmacophore queries. Lipinski rules of five were then applied to select only compounds with the best physiochemical properties for orally bioavailable drugs. 260 compounds were obtained and subjected to covalent docking-based virtual screening to determine their binding energy scores. The top three candidate compounds, which were shown to adapt similar binding modes as the reported covalent ligands were selected. The mechanism and stability of binding of these compounds were confirmed by 100 ns molecular dynamic simulation followed by MM/PBSA binding free energy calculation. The identified compounds can facilitate the rational design of novel covalent inhibitors of MERS-CoV 3CL(pro) enzyme as anti-MERS CoV drugs. Coronaviruses infect most vertebrate species, including humans. In humans, they generally cause mild-to-moderate respiratory tract infections [1] . However, viruses belonging to this group namely Severe Acute Respiratory Syndrome Coronavirus (SARS-CoV), MERS-CoV, and recently identified Severe Acute Respiratory Syndrome-related Coronavirus (SARS-CoV-2, which causes COVID-19 disease) have been shown to cause fatality in humans [2, 3] . MERS-CoV was discovered in 2012 in the Middle East [4, 5] . It is transmitted from camels to human and humanto-human transmission has been reported [6, 7] . As of January 2021, a total of 2566 confirmed cases of MERS, including 882 related deaths (fatality rate: 34.4%), have been reported. The majority of the cases were reported from Saudi Arabia (2121 cases), including 788 deaths (fatality rate: 37.1%) [8, 9] . MERS-CoV causes severe pneumonia in most cases of infection, quite similar to SARS. The patients are put on respiratory support and, unlike SARS, also monitored for renal failure [5, 10] . MERS-CoV belongs to the genus β-coronavirus (Subfamily: Coronavirinae; Family: Coronaviridae; Order: Nidovirales), and it is closely related to bat coronavirus [5, 11] . MERS-CoV is a singlestranded RNA genome of positive polarity and features the largest RNA genome known to date (~30 kb) [12] . The entry of MERS-CoV begins with the binding of its spike protein (S) on the dipeptidyl peptidase 4 (DPP4) receptor on the host cell surface [5, 11] . Upon entry of the virus particle, it is decoded and generate large polyproteins pp1a (4382 amino acids) and pp1ab (7073 amino acids) via translation of open reading frame (ORF 1a & 1b) [11, 12] . These polyproteins are then proteolytically cleaved by two viral proteases; 3-Chymotrypsin-like cysteine protease (3CL pro ) and papain-like protease (PL pro ), into at least 15 non-structural proteins (NSPs) [13] . 3CL pro cleaves the polyprotein at 11 distinct sites to generate many of the non-structural proteins, which are important in viral replication [13] . Thus, this protease plays a critical role in replication of virus [14] [15] [16] . The interruption of any of these replication processes would become a potential molecular target for antiviral drug development. 3CL pro shows weak dimerization and is believed to be dimerized in the presence of its substrate [17] . Structure-based activity studies and various high throughput studies have identified distinct inhibitors of SARS-CoV 3CL pro [16, [18] [19] [20] [21] [22] . Many of these inhibitors have been tested against MERS-CoV. However, there are no studies that show that SARS-CoV inhibitors have a distinct effect against MERS-CoV [23] . Thus, it is essential to identify novel inhibitors of MERS-CoV 3CL pro . Some peptidomimetic compounds could inhibit the dimerization of MERS-CoV 3CL pro [13, 24, 25] . The inhibitors reported to date comprise: N3 (IC 50 0.288 μmol/l, EC 50~0 .3 μM) [26] , 5-chloropyridyl esters GRL-001 [27] , benzotriazole derivatives [17] , dipeptidyl derivative GC376 ((EC 50 1.6 μM) [28] , pyrazolone derivatives (EC 50 5.8 μM to 7.5 μM) [25] , pyrrolidinone derivatives (EC 50 1.7 μM to 4.3 μM) [29] , pyrrolidinone-based peptide GC813 (EC 50 0.5 μM to 0.8 μM) and piperidine embodied peptidomimetics [24] . These inhibitors validated the scope of targeting this enzyme with covalent compounds with novel chemical scaffolds for strong anti-MERS-CoV activity. Moreover, 3CL pro is highly conserved across coronaviruses, therefore, there is a possibility for the identification of compounds that could have broad spectrum anti-viral activity [30] [31] [32] [33] . Covalent inhibitors contain electrophilic warheads including epoxide, aziridine, ester, ketone, α, β-unsaturated carbonyl, acetylene, nitrile, sulfur tethers, etc. that react and form covalent bond with nucleophilic residues like serine, cysteine etc. in the active site of the enzymes/proteins [32, 34] . They have strong target affinity and if the reactivity of the electrophilic warhead is controlled, these inhibitors can provide better therapeutic options due to their outstanding pharmacodynamic properties [32, 35] . The small molecule database could be utilized in either ligand-based or structure-based virtual screening for an effective identification of inhibitors [30] . In this contribution, a virtual screening approach based on a pharmacophore modeling followed by a covalent docking, molecular dynamic simulations and MM/PBSA binding free energy analyses were utilized to explore potential small covalent inhibitors of MERS-CoV 3CL pro enzyme as anti-MERS-CoV drugs. The workflow illustrating the methodology for the identification of novel MERS-CoV 3CL pro enzyme inhibitor is shown in Fig. 1 . 2.1. Establishment of a structure-based pharmacophore model LigandScout 4.3 program [36] was utilized to construct a pharmacophore model based on two X-ray structures of MERS-CoV 3CL pro in complex with covalent inhibitors. The PDB files of both structures, 5WKK and 5WKJ, with X-ray resolution of 1.55 and 2.05 Å, respectively, were obtained from protein data bank. Two independent pharmacophore hypotheses were initially generated from the 5WKK and 5WKJ of MERS-CoV 3CL pro structures in complex with GC813 and GC376 inhibitors, respectively, using a default setting [24] . The shared pharmacophore features along with exclusion volumes were then excreted to construct the final pharmacophore model that was used as 3D query for pharmacophore-based virtual screening of covalent chemical database. The developed 3D pharmacophore model was qualitatively validated using Receiver Operating Characteristic (ROC) method implemented in LigandScout software. The analysis was carried out by screening the model against two sets of compounds: 9 known actives and 512 decoys. The decoy set was generated from the active molecules using DUDE decoy online generator [37] . The virtual screening was performed with three covalent focused chemical databases consisting of a total 31,642 compounds. The libraries include Enamine covalent library (ww.enamine.net), ChemDiv covalent inhibitor library (http://www.chemdiv.com) and Life chemicals cysteine focused covalent inhibitor library (www.lifechemicals.com), composed of 21,969, 8293 and 1383 compounds, respectively. Discovery Studio Visualizer was used to combine all three databases in one SDF file [38] . The duplicates in the chemical database were removed by using OpenBabel (ver. 3.0.0) [39] . As a result, a covalent library of 28,790 compounds was prepared. A structure-based pharmacophore virtual screening was performed using LigandScout 4.3. The hits were ranked based on the pharmacophore fit scores and the only compounds that meet the allpharmacophore query were considered for further investigation. To reduce the number of candidate compounds, Lipinski rule of five [40] in LigandScout 4.3 program tool were applied. Compound with the following parameters is considered for docking-based virtual screening; molecule weight < 500 g/mol, LogP <5, Hydrogen bond donors <5, hydrogen bond acceptors <10 and rotatable bonds <10. To filter the candidate compounds further, the obtained 260 candidate compounds were subjected to covalent docking based virtual screening (CovDock-VS) against the structure of MERS-CoV 3CL pro (PDB: 5WKK). Molecular docking studies were carried out by using Maestro program, version-10.5 (GUI of Schrodinger 2017). X-ray structures of MERS-CoV 3CL pro having PDB ID of 5WKK with resolution of 1.55 Å was chosen for covalent docking [24, 41] . Protein Preparation Wizard of Maestro was used for the preparation and energy minimization of the crystal structure. Hydrogen atoms were added, water molecules were deleted, and the tautomeric and protonation states of Asp, Glu, Arg, Lys and His amino acids were adjusted and finally energy of the crystal structure was minimized by OPLS2005 force field. The ligands were prepared by LigPrep facility in Schrodinger. Various stereoisomers, tautomers, and the protonation and ionization states of ligands were generated at pH 7.4 using ionizer, as a result of which a library of 435 ligands were constructed. Finally, energy of the ligands was minimized by OPLS2005 force field. Before molecular docking, the reactive amino acid residue, Cys148 was mutated to alanine, to dock the prereactive forms of ligands. The ligands generated by LigPrep were docked into the receptor grid in the standard Glide XP mode. The reactive functional groups of the ligands were constrained within 5 Å of the C β atom of the reactive amino acid residue. Post-docking minimization was performed and up to three best poses per ligand were written out. All poses were examined manually, and the pose of the best docking score was retained unless noted otherwise. Then, the amino acid was mutated back to Cys148 and sampled using Prime VSGB2.0 with OPLS2005 force field. Covalent docking was performed using the CovDock application with Cys148 as the nucleophilic residue, which underwent a conjugate addition to an alkyne (carbonyl-activated) as preset by CovDock. The covalent bond was formed for ligand poses having reactive functional groups within 5 Å according to the reaction specified. The selection and ranking of the ligands were done on the basis of the Glide Scores of the binding mode of pre-reactive complexes. The available QikProp tool in Schrodinger was used for finding the drug-likeness properties. Based on Lipinski's rule of five, the properties that have been considered were molecular weight (MW), hydrogen bond donor (HBD), hydrogen bond acceptor (HBA), lipophilicity (log P), rotatable bonds and aqueous solubility (QP log S). The ADME-T (absorption, distribution, metabolism, and excretion -toxicity) profile is very essential for predicting the pharmacokinetics properties of compounds [42] . ADME-T properties of hit compounds were determined using QikProp module of Maestro program, version-10.5 (GUI of Schrodinger 2017). A 100 ns molecular dynamic (MD) simulation on covalently docked structures was performed using GROMACS 2018.1 package [43] with the OPLS-AA all-atom force fields [44] . MD simulations were carried out in the presence of ligands covalently bound to Cys148. The topology parameters of ligands were generated using SwissParam [45] and LigParGen [46] web servers. Each of docked-ligand and protein complex was then solvated in a triclinic box with TIP3P water molecules [47] with at least 1 nm distance from the docked-ligand and protein system. The system was further neutralized by adding appropriate numbers of counter ions. Periodic boundary conditions were applied during the MD simulation. A steepest decent algorithm of a maximum step size 0.01 nm with a tolerance of 239 kcal/mol/nm was used for energy minimization. A LINear Constraint Solver (LINCS) algorithm was used for bond lengths constrained. Electrostatic calculations were achieved using particle mesh Ewald method. After, energy minimization, the system was equilibration using canonical ensembles NVT followed by isothermal-isobaric ensemble (NPT) for 100 ps. All the simulations were carried out at constant temperature (300K) and pressure (1 atm). The production run was 100 ns, with trajectories generated every 2 fs and saved every 2 ps. All the preliminary analyses such as root mean square deviation (RMSD), root mean square fluctuations (RMSF), radius of gyration (Rg), and hydrogen bond analysis were carried out by GROMACS analysis programs. 2.6. MM/PBSA binding free energy calculation MM/PBSA binding free energy analysis was performed using G_mmpbsa module of GROMACS v.5.1.4 [48] . Briefly, Total 5000 snapshots from the last 20 ns stable simulation trajectories of every system were extracted with the interval of 2 ps for the calculation of MM/ PBSA. Following equations have been used to calculate binding free energy: Representations in the above equations are: ΔG bind = binding energy; -TΔS = entropic energy; ΔE MM = molecular mechanics potential energy [van der Waals (ΔE vdw ) energy + electrostatic (ΔE elec )]; ΔG sol = solvation free energy [polar solvation energy (ΔE polar ) + non-polar solvation energy (ΔE SASA ). MM/PBSA analysis has been widely used in binding free energy calculations for anti-viral inhibitors [32, [49] [50] [51] [52] [53] [54] [55] . MERS-CoV infection is recognised as a global health issue. Currently, there is no treatment or vaccine to combat or prevent MERS-CoV infection. MERS-CoV 3CL pro has emerged as an attractive target to inhibit replication of this virus [5, 56] . Crystal structures of this protease in complex with dipeptidyl and tetrapeptidyl peptidomimetic inhibitors, namely GC813 and GC376, respectively, showed that both compounds formed a tetrahedral hemithioacetal upon reaction with the catalytic Cys148 residue in the active site [24] . This shows the viability of targeting the active site Cys148 with covalent inhibitors. In fact, targeting protease enzymes with covalent inhibitors is an effective general strategy in antiviral drug discovery [32, 34] . In this context, the main objectives of this study were to use a pharmacoinformatic approach and molecular dynamic simulation to explore novel covalent inhibitors of MERS-CoV 3CL pro . This approach has been successfully used to discover new drugs in a time-and cost-effective manners [55, 57, 58] . Pharmacophore describes the three-dimensional arrangement of essential steric and electronic features for optimal binding of a ligand to a macromolecule. The design of pharmacophore model could be either structure-based or ligand-based depending on the available target or known ligands information, respectively [57, 59, 60] . In this study, the structure-based pharmacophore model was constructed from the shared features of two independent pharmacophore models that was designed from two X-ray structures of MERS-CoV 3CL pro in complex with two potent and covalent inhibitors namely GC813 and GC376. GC813 and GC376 inhibit the MERS-CoV 3CL pro enzyme at 0.5 and 0.9 μM, respectively ( Fig. 2A) [24] . Both compounds form a covalent tetrahedral adduct with the cysteine 148 (Cys148) of the catalytic dyad His41-Cys148 (Fig. 2B & C) . Overlapping both ligands suggested that they both adapt the same binding mode within the active site of MERS-CoV 3CL pro enzyme (Fig. 2D) . Using crystal structures of ligands in complex with their targets for generating the structure-based pharmacophore model may post the efficiency of the designed pharmacophore model [61] . The first pharmacophore model was designed from the interaction of GC813 with MERS-CoV 3CL pro enzyme. The pharmacophore was composed of one residue binding point which is the site for covalent binding group, one hydrophobic, three hydrogen bond donors and one hydrogen bond acceptor site (Fig. 3A) . The second pharmacophore model was generated from the binding of GC376 to MERS-CoV 3CL pro . This model was consisted of one residue binding point, three hydrophobic, five hydrogen bond donors and one hydrogen bond acceptor group (Fig. 3B ). Using the LigandScout 4.3 program, the shared features along with the exclusion volumes, which are essential features for the determination of the binding pocket overall shape, were extracted to result in the final pharmacophore model. This pharmacophore model was composed of one residue binding point, one hydrophobic, two hydrogen bond donors and one hydrogen bond acceptor group (Fig. 3C ). This pharmacophore model was then validated theoretically for its sensitivity and specificity to recognize true-active and -inactive ligands by screening a database composed of 9 known potent 3CL pro inhibitors and 512 decoy compounds using the receiver operating characteristic curve (ROC) method. A ROC curve illustrates the rate for the ability of a model to retrieve the true positive (actives) on Y-coordinate plotted against the rate for its ability to retrieve false positives (decoys) on Xcoordinate. An ideal curve would increase along the Y-axis until it reaches 1, which is the maximum true positive rate, then continues horizontally to the right shown that the hit list contains only the active molecules in the dataset. Two parameters of ROC curve were calculated at different fraction of the model-ordered database (1,5;10 and 100%) including the area under the curve (AUC), which shows the capability of the model to distinguish between true-active and decoy molecules as well as the enrichment factor (EF) which represents the number of true-active compounds found by using the generated pharmacophore mode [62] . The calculated ROC parameters for the generated model were as follows: AUC 1,5,10,100% 1.0, 1.0, 1.0, and 0.83, respectively; and EF 1;5;10;100%, 57.9; 15.8; 15.8; and 15.8, respectively (Fig. 3D) . These values complied with the requirements for AUC 1,5,10,100% all to be between 0.5 and 1 and EF 1;5;10;100% all to be >1 validating the ability of developed 3D model to distinguish the active molecules from decoys, therefore, its reliability to carry out the pharmacophore-based screening [63] . A focused covalent chemical library containing a total of 28,790 compounds was used for virtual screening against MERS-CoV 3CL pro enzyme. The compounds possess numerous reactive groups toward cysteine residue such as α,β-unsaturated ketones, α-chloracetamides, phenylsulphonate esters, vinylsulfonamides, acrylamides, acrylonitriles, aminomethyl methyl acrylathes, methyl vinylsulfones, epoxides, activated acetylenes and sulfonyl fluorides. This chemical database was initially used to carry out the pharmacophore-based virtual screening against the generated pharmacophore model. Among the screened compounds, 378 compounds were found to meet all the pharmacophore queries. The candidate compounds were ranked based on LigandScout pharmacophore fit score, which reflects to how much extent a compound fits the pharmacophore features. The pharmacophore fit scores for the obtained candidate compounds are ranged from 57.4419 to 55.0652. To reduce the number of candidate hits, Lipinski rule of five in LigandScout 4.3 program were applied. This step is necessary to evaluate the drug-likeness based on the physicochemical properties of the molecule. Among the 378 compounds, 260 candidate compounds were obtained and considered for the dockingbased virtual screening. To filter the candidate hits further, a covalent docking-based virtual screening was performed against the active site of MERS-CoV 3CL pro enzyme using CovDock-VS (Maestro ver. 11, Schrodinger). The 260 Table 1 The name, chemical structure, pharmacophore fit-score, Glide XP score, CovDock score and binding interaction of the hit compounds. Chemical structure (Chemical name) Pharmacophore fit-score (Table 1) . Interestingly, these compounds possess activated acetylene groups which are involved in a covalent interaction with the Cys148 of the catalytic dyad. The active site of MERS-CoV 3CL pro enzyme has a catalytic Cys148-His41 dyad and an extended binding site [24] . Before the molecular docking, the protocol was initially validated by independent redocking of the co-crystallized ligand, GC813 into the active binding site. The latter step was conducted to examine the ability of the docking protocol to produce the bioactive conformation. Accordingly, the docked pose with the lowest binding energy score adapted a binding mode similar to that of the co-crystallized ligand (Fig. 4A ). In addition, the distance between the reactive bisulfite group and the catalytic Cys148 residue was returned. The catalytic Cys148 was 1.8 Å from the bisulfite group of the co-crystallized ligand and 3.3 Å from the bisulfite group of the docked pose. These results validate the robustness of the docking protocol. The best binding poses for MA69, MA120 and MA152 are depicted in (Fig. 4B ). These candidate compounds adapted similar binding modes within the active site of MERS-CoV 3CL pro enzyme. The reactive acetylene groups within the chemical structures of these compounds formed covalent bond with the catalytic Cys148 ( Fig. 4C-E) . The reactive groups are oriented toward the catalytic Cys148-His41 dyad within the active site. Therefore, these compounds have high potential to inhibit the MERS-CoV 3CL pro enzyme. The selected hit compounds were then evaluated for the druglikeness and ADME-T properties by QikProp module of Schrodinger. From pervious filter, these hits follow the Lipinski's rule of five. These compounds have molecule weights (MW) < 500, log P < 5, hydrogen bond donors (HBD) <5, hydrogen bond acceptors (HBA) < 10. In addition, the number of rotatable bonds was <7, which is an essential parameter for good oral bioavailability. The value of QPlog S is an indicator of aqueous solubility of compounds with an acceptable range from −9.5 to 0.47. The lower the QP log S value, the higher the solubility and drug absorption profile. The identified hit compounds were found to perfectly obey the Lipinski's rule of five with ideal rotatable bonds numbers and QPlog S values, essential for orally active drugs [42] . Ghose et al. derived a guideline for the selection and optimization of CNS and non-CNS orally active compounds, through the analysis of 35 physicochemical and pharmacokinetic properties of 317 CNS and 626 non-CNS oral drugs [42] . Interestingly, the physicochemical and pharmacokinetic profile of compounds MA69, MA120, and MA152, was found to be within the limits as proposed and presented in the Table 2 . To confirm the screening and docking results, the complex of the three hit compounds with MERS-CoV 3CL pro were then subjected to 100 ns MD simulations to study their binding stability. The advantage of MD simulation is that the ligand-protein complex is studied under an environmental condition similar to that of human cell condition with respect to temperature, pressure, solvent and ions. Therefore, the data obtained from MD simulation could provide deep insight into the mechanism, dynamic and natural of ligand-protein interaction [65] . Over the course of 100 ns time, the ligand-protein complexes were analysed based on their root mean square deviation (RMSD), root mean square fluctuations (RMSF), radius of gyration (Rg), and hydrogen bond analysis. 3.5.1. Root mean square deviation (RMSD) and root mean square fluctuations (RMSF) To obtain the equilibrium time of each simulated protein-ligand complex during MD simulation, the backbone's root mean square deviation (RMSD) was calculated. RMSD plot is usually used to evaluate the time required for a system to reach structural equilibrium and to estimate the duration of running a simulation. The RMSD is a valuable parameter to estimate shifts or changes in molecular conformation. The Abbreviations: QL, qualifying lower limit; PL, preferred lower limit; QU, qualifying upper limit; PU, preferred upper limit [QL, PL, QU and PU values for non-CNS drug criteria were taken from [64] ]. RMSD values of simulated complexes including the reference were suddenly raised due to the sudden change in the structure condition, which is related to the protein crystallization method. The latter effect was expected, as in the crystal structure the protein is rigid and when it gets solvated in water box it returns back its dynamic motion [66] . The backbone RMSD values of the simulated systems revealed that an equilibration was achieved after~20 ns for all complexes with respect to the reference; 3CL pro -GC813 complex (Fig. 5A ). Unlike other complexes, there was an observed increase in the RMSD of compound MA69-3CLpro complex up to 0.45 Å at 84 ns (Fig. 5A ). This could be a result of adapting a new conformation within the active site by the ligand. The averaged RMSD values of GC813, MA69, MA120 and MA152 complexes, for the last 80 ns were 0.21 ± 0.02-, 0.21 ± 0.03-, 0.20 ± 0.02and 0.17 ± 0.02 nm, respectively. The RMSF was calculated to assess the contribution of each residue to the complex fluctuation. The RMSF of each residue within the MERS-CoV 3CL pro structure were calculated. Residues with high fluctuation pattern such as Glu197 and Glu277 are located in the loop regions away from the active-catalytic site. The RMSF values suggested normal backbone fluctuations behaviour indicating stable hits binding (Fig. 5B ). The radius of gyration (Rg) which is an indicator of the protein compactness suggested normal behaviour in the protein structure with average values of 2.18 ± 0.01-, 2.20 ± 0.01-, 2.19 ± 0.01-and 2.18 ± 0.01 nm for GC813, MA69, MA120 and MA152 complexes, respectively (Fig. 6A) . The hydrogen bond is a major factor responsible for stable conformation's maintenance of the protein and ligands [67] . H-bond analysis was performed to examine the H-bonds for ligands complexes over the 100 ns of the simulation time and presented in Fig. 6B . The Ligands' complexes showed up to two MA69, up to four MA120 and up to three MA152, during 100 ns of the simulation with respect to GC813 which showed up to 7H-bonds. The Molecular Mechanics Poisson-Boltzmann Surface Area (MM/ PBSA) is an effective and dependable approach to calculate binding free energies of small inhibitors to their protein targets [32, 68] . The binding free energies and energy components of all 3 systems are listed in Table 3 . The calculated bind free energy of complex MA69 is −27.75 kcal/mol, MA120 is −19.33 kcal/mol and of complex MA152 is −25.21 kcal/mol. For all the 3 complexes, ΔE vdw and ΔE elec were favourable and showed stable interaction between ligands-protein. The 100 ns molecular dynamic simulation study reveals that the binding of these compounds to the MERS-CoV 3CL pro were stable throughout the simulation course. Collectively, the identified compounds may have high potential of inhibition of this critical enzyme for replication of MERS-CoV. In this study, we identified novel irreversible inhibitors of MERS-CoV 3CL pro as potential anti-MERS-CoV drugs by integrated computational approaches including pharmacophore modeling, covalent docking, molecular dynamics simulation and binding free energy calculation analyses. These compounds carrying electrophilic acetylenes groups as reactive binding points interact covalently with the key catalytic Cys148 of the MERS-CoV 3CL pro with high binding affinity forming stable acetylene adducts. These compounds with good drug-like properties, may serve as seeds to rational development of potent irreversible MERS-CoV 3CL pro . 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. of the International Anti-Cancer Drugs The authors would like to thank Prince Sattam Bin Abdulaziz University, AlKharj, Saudi Arabia for providing necessary facilities to carry out this research.