ABBV-744 and Onalespid as Potential Inhibitors of SARS-CoV-2 main Protease Enzyme: A Promising Therapeutics in COVID-19?

A new pathogen severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has spread worldwide and become pandemic with thousands new deaths and infected cases globally. To address coronavirus disease (COVID-19), currently no effective drug or vaccine is available. This necessity motivated us to explore potential lead compounds by considering drug repurposing approach targeting main protease (M pro ) enzyme of SARS-CoV-2. This enzyme considered to be an attractive drug target as it contributes signi�cantly in mediating viral replication and transcription. Herein, comprehensive computational investigations were performed to identify potential inhibitors of SARS-CoV-2 M pro enzyme. The structure-based pharmacophore modeling was developed based on the co-crystallized structure of the enzyme with its biological active inhibitor. The generated hypotheses were applied for virtual screening based PhaseScore. Docking based virtual screening work-ow was used to generate hit compounds using HTVS, SP and XP based Glide GScore. The pharmacological and physicochemical properties of the best hit compounds were characterized using ADMET. Molecular dynamics simulations were performed to explore the binding a�nities of the considered compounds. Binding studies revealed that compound ABBV-744 binds to the M pro with strong a�nity ( (cid:0) G bind -45.43 kcal/mol), and the complex is more stable in comparison with other protein-ligand complexes. Our study classi�ed three best compounds which could be considered as promising inhibitors against main protease SARS-CoV-2 virus.


Introduction
Coronavirus disease 2019 (COVID-2019) outbreak is a global pandemic caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) which initially diagnosed in Chinese patients of Hubei's Wuhan city in early December 2019 [1].SARS-CoV-2 disclosed a close genetic resemblance to the severe acute respiratory syndrome coronavirus (SARS-CoV) that already triggered an epidemic in 2003 [2].COVID-19 has been declared a global health disaster by World Health Organization (WHO) on 30 January 2020 as the disease hastily transmitted human-to-human and affected more than 170 countries across the world [3].The existing condition is extremely increasing; therefore, the overall asperity of this disease persist to be serious.The infection rate of SARS-CoV-2 is higher (10-12%) in comparison with its mortality rate (5.4%) [4].The most distinctive indications of COVID-19 patients are high fever, cough and excessive respiratory sickness that required urgent intensive care facility.Currently, there is no applicable and precise medication for the treatment of COVID-19, however, many drugs and vaccines are under clinical trials.The only practical approach available is the repurposing of existing antiviral drugs as these drugs have already been tested for their toxicity [5].Still there is a prompt requirement to make substantial efforts to advance therapeutic interventions against CoV infections.
CoVs are single-stranded positive-sense RNA viruses belongs to the family of Coronaviridae.These viruses can be categorized into four species: alpha, beta, gamma and delta.The recent SARS-CoV-2 is from beta genus and is usually identi ed to affect commonly humans [6].The RNA genome length of this virus is about 27-32 Kb encoding both structural and non-structural proteins.The structural proteins are membrane (M), envelope (E), nucleocapsid (N), hemagglutinin-esterase (HE) and spike (S) proteins contribute notably in viral transmission and replication in the host cells [7].Among which 3C-like protease (3CLpro) protease plays critical role in virus replication and transcription, thus studied as potential drug targets.Providing informative knowledge regarding the enzyme inhibition will be valuable in tailoring effective and selective inhibitors of 3CLpro that will ultimately be associated to report COVID-19 as this protease is imperative for virus assembly and reproduction [8].
The main protease (M pro ) is a quintessential enzyme which contributes signi cantly in the life cycle of SARS-CoV-2 and inhibition of M pro enzyme activity would block viral replication.Since no human proteases with a similar speci ed cleavage are characterized, thus the potential inhibitors are unlikely to be toxic.The M pro enzyme consists of an asymmetric unit including 306 amino acid residues with CYS145 and HIS41 catalytic dyad in the active site [9][10] (Figure 1).
Experimental observations by Zhang et al. [10] demonstrated the half maximal inhibitory concentration (IC50) value 0.67 mM for an a-ketoamide inhibitor (ak-13b) as a potent antiviral inhibitor against M pro enzyme.Zhang et al. proposed the inhibition mechanism through the nucleophilic attack of the catalytic CYS145 to the α-keto group of the inhibitor forming thiohemiketal.This thiohemiketal group is stabilized by formation of several hydrogen bonds with the active residues of HIS41, GLY143, SER144 and CYS145, Figure 1.This experimental observation proposed that α-ketoamide inhibitors interacted with the catalytic center of proteases through two hydrogen bond interactions whereas aldehydes and Michael acceptors only formed one hydrogen bond into the catalytic center of the target proteases [11] [12].
According to the aforementioned experimental structural information, we have chosen SARS-CoV-2 M pro as a target enzyme to accelerate the prompt hunt of antiviral drug repurposing with clinical potential to gain a short-term and speci c solution to treat COVID-19 patients.
To address this challenge, an speci c library of anti and pro-viral agents including FDA approved drugs, compounds in clinical trials and preclinical compounds having inhibitory activity between 10-100 nM range against SARS-CoV-2 was considered for drug repurposing to attain immediate and precise results [13].In this study, we developed an integrated approach of drug discovery integrating 3D structure-based pharmacophore modelling, virtual screening of the considered library, molecular docking work ow, ADMET pharmacological analysis and molecular dynamics (MD) simulations to nd potential antiviralagents for the therapeutic management of COVID-19.This scheme will provide an informative insight into the exploration of potent antiviral drugs which, could help in progressive attempts in the therapeutics of COVID-19.

System preparation
The 1.95 Å crystal structure of SARS-CoV-2 main protease (M pro ) in complex with inhibitor was extracted from the Protein Data Bank (PDB ID: 6Y2F) [10].The structure of the enzyme was pre-processed, minimized and re ned using the Protein Preparation Wizard [14] implemented in Schrödinger suite.This involved eliminating crystallographic waters, missing hydrogens/side chain atoms were added, and the appropriate charge and protonation state was assigned to the receptor structure corresponding to pH 7.0 considering the appropriate ionization states for the acidic as well as basic amino acid residues.
Structure was subjected to energy minimization using the OPLS-2005 force-eld [15][16] with a root mean square deviation (RMSD) cut-off value of 0.30 Å to relieve the steric clashes among the residues due to the addition of hydrogen atoms.
The preparation of the crystalized inhibitor, ak-13b and the candidate compounds were achieved using LigPrep module of Schrodinger Suite which undertakes hydrogens atom addition, amending realistic bond lengths and angles, accurate chiralities, ionization states, tautomers, stereo chemistries, and ring conformations.Partial charges were assigned to the structures using the OPLS-2005 force-eld [15] and the subsequent structures were imperiled to energy minimization until their average RMSD reached 0.001 Å.The ionization state was set at the neutral pH using Epik ionization tool [17].

Preparation of Inhibitor-like Ligand Library
We have retrieved the 75 candidate compounds by Gordon et al. [13] as these compounds were selected based on their experimental anti-viral activity against SARS-CoV-2 of coronavirus.The total 75 compounds were considered for further virtual screening analysis.The list of the compounds has been displayed in Table S1.

Identi cation of 3D-Pharmacophore Hypotheses
For the structure-based pharmacophore modeling Schrodinger PHASE module [18] was used with the default set of seven chemical features-hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic contacts (H), negative ionizable (N), positive ionizable (P), and aromatic ring (R) to create the utmost illustrative features of the M pro active sites.The seven 3D-features were generated using Hypothesis Generation for Energy-Optimized Structure Based Pharmacophores considering the omitted volumes within 5Å of the re ned ligand for the enzyme [19].Pharmacophore features were chosen based on the crucial interactions with the key residues of the enzyme accommodated the inhibitor.The resulted pharmacophore features comprise the functional groups included in their bioactivity of targeted enzyme.The omitted volumes contain all atoms within 5Å of the re ned ligand for the target.

Screening of M pro Inhibitors
All acquired seven 3D-pharmacophore features were exported and used for PHASE-based virtual screening to screen the candidate library of 75 compounds which were retrieved from Gordon et al [13] recent experimental work.Out of 75 compounds, 43 were generated based on the highest PHASE screen score and matched ligand sites (Table S2).Both the quantity and quality of feature matching is taken into account in the Phase-Screen-Score factor.

Docking-based virtual screening
Molecular-docking-based virtual screening was performed using Glide-Based virtual screening work ow of Maestro 11.6 to prioritize the hit candidate compounds that strongly bind to M pro enzyme [20].
Receptor grid was created as center coordinates (X = 9.81 Y = -1.47Z = 20.51)using two cubic boxes having a mutual centroid to systematize the calculations: a larger enclosing and a smaller binding box with dimensions of 24 × 24 × 24 Å and 18 × 18 × 18 Å, correspondingly.The grid box was centered on the centroid of the ligands in the complex, which was adequately large enough to search a superior region of the enzyme structure.All the chosen ligands were docked by using a three docking protocols of Glide [20] which starts with "High throughput Virtual Screening" (HTVS) followed by "Standard Precision" (SP) and then by "Extra-Precision" mode (XP).Finally, 43 input compounds were assessed using Docking-Based Virtual Screening and ltered to nal three optimized lead compounds based on XP-GScores.

ADMET Properties Assessment
Schrodinger QikProp 5.6 module was used to calculate absorption, distribution, metabolism, excretion and toxicity (ADMET) properties of the considered compounds to produce the ADMET associated descriptors.This protocol predicts noteworthy physicochemical and pharmacokinetic-based descriptors of the compounds based on Lipinski's rule of ve [21] [22].ADMET properties of the top three compounds and one control inhibitor were considered and studied using QikProp 5.6 module and the best three compounds were chosen for nal molecular dynamics (MD) simulations.

MD Simulations
MD simulation considered to be the most essential approach in understanding the fundamental structure and function of biological macromolecules.This method helps in nding the underlying dynamics and how it is connected to enzyme's biomolecular function [23][24] [25].AMBER 18 [26] simulation package was used to execute 200 ns MD simulations on all the prepared complexes using (Graphics Processing Unit) GPU accelerated version of Partial Mesh Ewald Molecular Dynamics (PMEMD) simulations [27].The ff99SB [28] and the general AMBER force elds (GAFF) [29][30] were employed to parametrize the enzyme and the considered ligands using LEaP implemented in Amber 18.
The ANTECHAMBER module was used to assign atomic partial charges for the ligands employed in General Amber Force-Field (GAFF).The system was solvated using the TIP3P [31] explicit water in a cubic box with 8 Å box edge.The Na + counter ions were added to randomly to neutralize the complex.The partial Mesh Ewald (PME) [32] method was used to account the long-range electrostatic forces using cutoff of 12 Å, and the SHAKE algorithm [33] was used to constrain all the hydrogen atoms bonds.
Energy minimizations were performed in two stages with 2500 steps of steepest decent minimization followed by 2500 of conjugated gradient to remove the bad contacts.The rst stage was followed with a harmonic restraint of 500 kcalmol -1 A -2 on the solute molecule whereas, ions and water molecules were relaxed.On the second stage of minimization the restraints were removed and the whole system was relaxed.Each minimized complex was then gradually heated up from 0 K to 300 K for 200 ps to keep the solute using a weak harmonic restraint of 10 kcalmol -1 A -2 .The 50 ps density equilibration with weak restraints followed by the 500 ps constant pressure equilibration at 300 K were performed at constant pressure using Berendsen barostat [34].Ultimately, the production phase of 200 ns MD simulation was performed on all the complexes at a constant temperature of 300 K and constant pressure at 1 atm.

Post-dynamic trajectories analyses
The 200 ns MD trajectories were analyzed to calculate the RMSD of C α atoms, root mean square uctuation (RMSF) of each residue in the complex, radius of gyration (R g ), solvent accessible surface area (SASA), intramolecular and intermolecular hydrogen bond interactions using CPPTRAJ module [35] implemented in AMBER 18.Molecular visualizations and plotting were conducted using Maestro 11.6 and OriginPro 2018 software [36].

Principal component analysis (PCA)
PCA as an important tool for identifying the conformational changes of proteins was carried to describe the residual motions upon inhibitor binding of biomolecular complex [37].PCA generates highly correlated and anti-correlated uctuations derived from MD trajectories by applying dimensional reduction [38][39].The collective motions were studied using the positional covariance matrix C constructed based on the atomic coordinates and their corresponding eigenvectors.The eigenvalues and eigenvectors are de ned as the extent and the direction of motions, respectively [40][41].By the following equation, the matric elements of the positional covariance matrix C were determined: (see Equation 1 in the Supplementary Files) where q i and q j are the cartesian coordinates for the i, jth of Cα atom, and N is the number of Cα atoms.
To remove all translational and rotational movements, the average is calculated after superimposition with a reference structure using a least-square t procedure to excerpt the important motion from MD trajectories [42][43] [44].To derive the eigenvalues and eigenvectors, the symmetric matrix C is transformed into a diagonal matrix Λ of eigenvalues by an orthogonal coordinate transformation matrix T: (see Equation 2in the Supplementary Files) in which the eigenvectors correspond to the direction of motions relative to and each eigenvector associate with an eigenvalue that represents the total mean-square uctuation of the system along the corresponding eigenvector.CPPTRAJ module from the Amber 18 suite was used to perform the PC analysis and the porcupine plot of protein collective motions was created by NMWiz implemented inVMD [45].

Binding free energy calculations
The relative binding free energies were computed using Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) binding free energy method [46].Water molecules and counter ions were stripped using the CPPTRAJ module.The binding free energies (DG bind ) were calculated with the MM-GBSA method for each complex as below: (see Equations 3-8 in the Supplementary Files) The gas phase energy (DE gas ) is the sum of the internal (E int ), van der Waals (E vdW ) and Coulombic (E elec ) energies, (Eq.6).The solvation free energy is the pattern of polar (G GB ) and nonpolar (Gnonpolar, solvation) contributions (Eq.7).The polar solvation G GB contribution is estimated using the Generalized Born (GB) solvation model with the dielectric constant 1 for solute and 80.0 for the solvent.Conversely, the nonpolar free energy contribution was assessed using Eq. 8, where the surface tension proportionality constant, g, and the free energy of nonpolar solvation of a point solute, b, were set to 0.00542 kcal mol -1 Å -2 and 0 kcal mol -1 , respectively [47].The SASA is calculated by the linear combination of pairwise overlap (LCPO) model [48].

Selection of compounds
An speci c drug repurposing library of 75 anti and pro-viral agents including FDA approved drugs, clinical trials compounds and preclinical compounds having inhibitory activity between 10-100 nM range [13] against SARS-CoV-2 was considered as the input database for this in-silico study, Table S1.

Structure-based pharmacophore modeling
The comprehensive and accurate information of ligand interacting features can be obtained from structure-based pharmacophores based on three-dimensional structure of a target protein [49].The most common descriptors in pharmacophore modeling are H-bond donors, H-bond acceptors, positive and negative ionizable groups, lipophilic regions and aromatic rings.The most effective 3D structure-based epharmacophores were produced using the receptor-ligand pharmacophore generation protocol implemented in PHASE, which was executed for a co-crystal ligand inside the active pocket in order to determine possibly critical amino acids that are involved in ligand binding (Figure 2A).The generated e-pharmacophore for the considered enzyme showed seven main 3D-features including, H-bond acceptor, H-bond donor and aromatic rings.In each pharmacophore model, the red arrows represent hydrogen bond acceptor, blue arrow represents hydrogen bond donor and orange spheres represent an aromatic ring in Figure 2B.Numerous excluded volumes were also produced in the models to demonstrate the space balancing.The 3D pharmacophore features and 2D-chemical structure of 3b are presented in Figure 2B and 2C showing seven 3-D pharmacophore features as three donor hydrogen bonds, three acceptor hydrogen bonds and one aromatic ring sphere.

Virtual Screening of the candidate compounds
The obtained structure-based pharmacophore hypotheses of ak-3b-inhibitor in complex with M pro were used to screen the 75 candidate anti-viral compounds retrieved from recent experimental work by Gordon et al [13] (Table S1).These compounds were screened based on PHASE screen score, matched ligand sites indices.A total of 43 compounds passed this lter and subsequently 43 compounds were ltered based on the created pharmacophore hypothesis.Molecules which have satis ed all the features of the pharmacophore model were considered as potential hits.The virtual screened compounds are presented in Table S2.

Docking-based virtual screening analysis
A total of 43 screened compounds obtained from virtual screening were considered for docking analysis using Glide module [50] of Schrödinger package.Three step wise ltering protocol were used for docking using HTVS where a total of 23 compounds (Table S3) were obtained followed by Glide SP where a total of 12 hits were generated (Table S4).Finally, the best ve compounds were ltered using Glide XP lead optimization docking while only one pose per ligand was retained (Table S5).The Glide GScore and the interacting binding residues of the ve lead compounds presented in Table 1.The ak-13b-inhibitor as well as the best optimized lead generated from XP docking were selected to map their potential interactions within the active pocket of SARS-CoV-2 M pro enzyme using molecular docking approach.This approach aids in understanding the optimized orientation of a ligand and its target protein by minimizing inclusive energies of the corresponding complexes.The estimated docking binding energy values of all three compounds and ak-13b-inhibitor are shown in Figure 3 and 4 and Table 1.
The Daunorubicin-M pro , Onalespib-M pro and ABBV-744-M pro docked complexes presented noteworthy binding a nities with the energy values of -9.33, -8.21 and -7.79 kcal mol -1 , respectively (Table 1).These three compounds contributed into the binding site interactions through the non-covalent interactions including hydrogen bond interactions, p -p interaction and p-Sulphur interactions present in the binding site of M pro enzyme (Figure 4).
The daunorubicin created one hydrogen bond and one p -alkyl interaction with the catalytic dyad CYS145 and HIS41, respectively.The other ve hydrogen bonds were formed by the hydroxyl (-OH) groups of daunorubicin with HIS164, Arg188, Thr190, ASP187 and Gln192 as presented in Figure 4A.ABBV-744 formed the interaction network of six hydrogen bonds with Arg188, Thr190, GLU166, HIS163, GLN189 and Gln192 with the hydroxyl group (-OH) of the compound, Figure 4B.The catalytic dyad CYS145, HYS41 and MET165 formed p-p stacking interaction and p-Sulphur interactions with ABBV-744.

ADMET analysis
Pharmacokinetic and toxicity features were predicted using QikProp module of Schrodinger for Daunorubicin, Onalespid, ABBV-744 and ak-13b-inhibitor.Outcomes of pharmacokinetic and toxicity study are illustrated in Table 2.The selected properties of the compounds are representatives of in uence metabolism, cell permeation, bioavailability and toxicity.The predicted central nervous system activity (CNS) of Daunorubicin, ABBV-744 and ak-13b inhibitor depicted as inactive whereas, Onalespib was presented as an active compound.The predicted human binding serum albumin (QPlogKhsa) of all compounds showed in the acceptable range.The estimated total solvent accessible surface area (SASA) of all three compounds and ak-13b inhibitor met the acceptable range: 300-1000.Predicted octanol/water partition coe cient (QPlogPo/w) showed in the acceptable range from −2 to 6.5 for all the ligands.The predicted aqueous solubility (QPlogS) for Daunorubicin, Onalespib and ak-13b inhibitor were in the acceptable range of −6.5-0.5 whereas ABBV-744 showed slightly low value.The predicted brain/blood partition coe cient (QPlogBB) for all these compounds showed in the acceptable ranges.The percentage human oral absorption for all the compounds met in the recommended range.Number of violations of Lipinski's rule of ve for all the compounds satis ed this rule for all the studied ligands.

Post-dynamics trajectories analysis
The structural variations within the enzymes structure is correlated with their biological activities.Any alterations or interference on enzymes structural integrity might have a substantial impact on its activity [39].The binding of inhibitors in uence the mode of action of enzymes that are comprised in disease pathways, thus there is a requirement to estimate the structural dynamics and conformational changes associated with the inhibitory activity of these inhibitors [51].
In this section, 200 ns MD trajectories regarding the four complexes, namely, Daunorubicin-M pro , Onalespib-M pro , ABBV-744-Mpro and ak-13b inhibitor-M pro as control model were analyzed.Different metrics and analysis were applied to investigate the stability and exibility of the complexes as well as the contribution of the studied compounds upon binding in terms of binding free energies.The 2D chemical structure of all the ligands considered for MD simulations have been displayed in (Figure 1 and Scheme 1).
Scheme 1: 2D chemical structure of top three hit compounds.(see Supplementary Files) The computation of a time variable with reference to an RMSD of C α atoms from generated trajectories was accomplished to investigate the consistency and e ciency of M pro in complex with ak-13b inhibitor and along with the top three compounds, Figure 1 and scheme 1.
The perturbations in the RMSD values as denoted in plot (Figure 5A) throughout the simulation time disclosed the possible conformational deviances in the enzyme structure upon ligand binding.As Figure 5A revealed, all the complexes were stabilized and attained convergence after almost 50 ns of simulation run.ABBV-744-M pro unveiled the lowest average RMSD of 2.45 Å, while Onalespid-M pro and Daunorubicin-M pro revealed average RMSD of 2.73 Å and 2.76 Å respectively.The ak-13b-inhibitor-Mpro unveils a perturbation of 2.85 Å as indicated in the plot.This evaluation proposed that any further analyses performed on the produced trajectories of all complexes were reliable.The RMSD plots indicated that ABBV-744-Mpro, Onalespid-M pro and Daunorubicin-M pro complexes exhibit the lowest deviation of C α backbone atoms suggesting that stability of these three compounds in icted higher constancy on M pro enzyme as compared to the ak-13b inhibitor-M pro complex.
To provide detailed insight into the structural uctuation and exibility of different regions of the amino acid residues of M pro enzyme upon binding of the selected compounds, RMSF values for Cα atoms were calculated from trajectories generated over 200 ns of MD trajectories.Hence, ligand binding to the enzyme could be investigated in relation to the modi cation in exibility in terms of RMSF values [52].To discover the stringency and elasticity in M pro residues upon binding of chosen compounds, RMSF values for C α atoms were estimated from trajectories produced from 200 ns of MD simulations run.As shown in Figure 5B, ABBV-744-M pro complex showed the least uctuations in the amino acid residues with 7.56 Å.
An average RMSF of 13.76 Å and 15.09 Å was spotted in complex Daunorubicin-M pro and Onalespid-M pro , correspondingly.The complex, ak-13b inhibitor-M pro disclosed an average of 14.11 Å that is remarkably greater than ABBV-744-M pro complex, signifying enhanced binding in comparison to the ak-13b inhibitor-M pro complex.This noteworthy decline might be coupled with structural inactivation that evidently con rmed as a result of signi cant binding of ABBV-744 compound in the active pocket of M pro enzyme.The reduced uctuation of amino acid residues might have favored M pro enzyme inhibition through compound ABBV-744.
The radius of gyration (Rg) parameter was assessed as the structural compactness index and its folding and unfolding behavior through the overall conformational variations in enzyme structure upon inhibitor binding.The mediocre values of R g for ABBV-744-M pro , Onalespid-M pro and Daunorubicin-M pro complexes were noted to be 41.55 Å, 41.87 Å and 42.08 Å, respectively.Figure 5C plots disclosed extremely minor changes in the compactness of the three compounds.The compound ABBV-744 exhibited a lowest R g in comparison with other two complexes and also with the control ak-13b inhibitor-M pro complex (42.36 Å), thus proposing improved compactness and enhanced binding with the M pro enzyme.All these patterns of conformational analysis are suggesting an improved stability, exibility and compactness of compound ABBV-744 in complex the M pro enzyme.
Solvent Access Surface Area (SASA) analysis was performed to de ne the activity of hydrophobic and hydrophilic amino acid residues and forces exposed to the solvent over 200 ns MD trajectories.The constant and accurate scheming of SASA is highly useful in the energetic evaluation of biological macromolecules [53].The interfaces among the hydrophobic native contacts inside enzyme structure is a noteworthy intermolecular interaction that effect enzyme inhibition.Hydrophobic interaction constructed between the non-polar residues corroborate the stability of the enzyme structure in solution by sheltering the non-polar residues inside the hydrophobic core distant from an aqueous solution [54].As displayed in Figure 5D, standard SASA values for all selected compounds have been measured during 200 ns MD trajectories.Average value of SASA for the compound ABBV-744-Mpro complex is 14230 Å 2 which was showed to the solvent system.Overall SASA values of 14303 Å 2 and 14426 Å 2 were prominent by Onalespid-M pro and Daunorubicin-M pro complexes, individually.The differences in SASA values for all the complexes during the simulation period corresponds with the folding and unfolding of enzyme structure.
The overall SASA values in the control complex was 14001 Å 2 , slightly less than ABBV-744-M pro complex.
The SASA assessment perceived in compound ABBV-744 bound complex additionally validated that ABBV-744 compound has better exposure to solvent and consequently favored the improved inhibitory activity of compound ABBV-744 over other complexes.

Hydrogen bond analysis
For overall conformation and stability of enzyme structure, we have measured the intramolecular and intermolecular hydrogen bond analysis (Figure 6).This analysis gives extreme understanding into binding mechanism of enzyme-ligand with detailed consideration [55].An average number of intramolecular hydrogen bonds in ABBV-744-M pro complex was noted to be 136 as displayed in Figure 6A.In compound Onalespid and Daunorubicin, the intramolecular hydrogen bonds were observed to be 139 and 140, respectively.
The number of intermolecular hydrogen bonds produced in the catalytic site of M pro enzyme notable to be 9-10 in ABBV-744 bound M pro complex.However, number of these bonds are more in Daunorubicin-M pro complex with 11-12 hydrogen bonds and less in Onalespid and ak-13b-inhibitor bound M pro with 7-8 hydrogen bonds as presented in Figure 6B.

Principal Component analysis
To qualitatively probe the impact of inhibitor's binding on the predominant conformational motion of each residue, [56][57] the concerted motions in Daunorubicin-M pro , ABBV-744-M pro and Onalespid-M pro and ak-13b-inhibitor-M pro complexes were studied using PC analysis based on the eigenvector.The scatter plots in It is evident in Figure 7A, the Onalespid, Daunorubicin and ak-13b-inhibitor with the trace covariance matric of 37.45 Å 2 , 37.67 Å 2 and 37.69 Å 2 , imposed highly uctuated anti-correlated effects as the negative values of 2D-scatter points into the protein, Figure 7B.Interestingly, in the case of ABBV-744 with the trace covariance matric of 37.64 Å 2 , the prominent correlated motions were observed with the least uctuations of the system upon ligand binding, Figure 7B.Thus, from the above observations Figure 7B, it was concluded that ABBV-744 induced least uctuations into the binding site upon binding than the variants complexes.

Mechanistic insights into binding a nity
To understand the impact of inhibitors upon complexation in terms of their binding a nities, MM-GBSA binding free energy method were utilized to calculate the binding free energies and their components of the complexes, Table 3.
Table 3. Binding free energies and its components for the three hit compounds:M pro and ak-ab-inhibitor:M pro using MM-GBSA method.The energy components are in kcal mol - As it is evident in Table 3, the total binding free energies (ΔG bind ) of Daunorubicin-M pro , Onalespib-M pro , ABBV-744-M pro and ak-3binhibitor-M pro were -36.65 kcal/mol, -37.13 kcal mol -1 , -45.43 kcal mol -1 and -20.92 kcal mol -1 , respectively.Accordingly, among all the studied complexes, ABBV-744-M pro and Onalespib-M pro depicted the most favorable of ΔG bind with lowest values of -45.43 kcal mol -1 and -37.13 kcal mol -1 .At this point, it is interesting to address the key contributions that each binding component can impose to the total binding free energies.
It is evident that amongst the studied complexes, the ΔG gas as the favorable contributing index into the total ΔG bind has the lowest values for ABBV-744-M pro (-98.65 kcal mol -1 ) and Onalespib-M pro (-62.10 kcal mol -1 ) complexes.This observation implies the most favorable contribution values of ΔE vdw and ΔE elec for ABBV-744-M pro (-54.45 kcal mol -1 and -35.20 kcal mol -1 ) and Onalespib-M pro (-43.55 kcal mol -1 and -18.55 kcal mol -1 ) into the total binding free energies.
Interestingly, this observation revealed another key contributing component of DG nonpolar for the complexes of ABBV-744-M pro (-5.81 kcal mol -1 ) and Onalespib-M pro (-4.45 kcal mol -1 ) which leads to lowest DG bind values for both complexes.According to the obtained energy components results, it could be inferred that the DE vdw , DE elec and DG nonpolar have the dominant contribution into the binding a nities for the selected complexes.

Per-residue decomposition energy analysis
The binding-free energy decomposition offers an immense understanding in account of enzyme-ligand complexes produced from the trajectories by MD simulations.To achieve this, we fragmented the total binding energies of complexes into each-residual involvement by per amino acid residue existing in the catalytic site of M pro enzyme to provide comprehensive identi cation of key contributing residues upon ligand binding as depicted in Figure 8.The interactions between catalytic site electro-negative and electro-positive residues develops ligand binding and its stabilization at the target enzyme.This creates an improved intermolecular binding that surges the binding a nity of the ligand in the active site.Active site residue Met165 in ABBV-744 contributed with the lowest ΔG bind with -4.22 kcal mol -1 however this residue has contributed with notably less ΔG bind of -2.50 kcal mol -1 and -2.45 kcal mol -1 in Onalespid and Daunorubicin bound complexes.The ΔG bind of another participating residues Gln189 was also observed to be lowest in ABBV-744 bound M pro complex with -3.94 kcal mol -1 however, it is slightly less in Onalespid-M pro complex with value of -3.45 kcal mol -1 and -2.48kcal mol -1 in Daunorubicin-M pro complex.Gln192 signi cantly contributed in the binding of ABBV-744 compound with ΔG bind value of -3.00 kcal mol -1 and observed a very minor difference with -2.99 kcal mol -1 ΔG bind in Onalespid bound M pro but showed a lesser ΔG bind of -1.50 kcal mol -1 in Daunorubicin compound.

4.
The necessity to control alarming COVID-19 pandemic made us to rationalize potential lead compounds that could be considered in clinical trials.Despite major investigations in the design and development of speci c drugs or vaccines, not much proven to be effective against COVID-19.This challenge motivated us to explore the drug designing approaches that could serve informative to combat this disease.In this report, we have performed 3Dstructure-based pharmacophore modeling followed by virtual screeningbased 3D-pharmacophore hypotheses of 75 compounds as potential antiviral agents retrieved from PubChem.Molecular docking work ow using HTVS, SP and XP protocols were used to generate the best hits and their corresponding docked poses.The Six best compounds generated based on their lowest docking binding a nities using XP were considered for ADMET prediction-based physicochemical and pharmacokinetic descriptors and MD simulations analysis.MD simulations approach revealed the two highly selective compounds namely, ABBV-744 and Onalespid possessed signi cant binding a nity and presumably inhibition of SARS-CoV-2 M pro enzyme.Based on our overall observations, compounds ABBV-744 and Onalespid could be recommended as potential lead for the therapeutic of COVID-19 patients.

Declarations Figures
Figure, essentially give a two-dimensional representation of the conformational changes occupied by the system.The gradual migration of the points in the PC1−PC2 scatter plots obtained using construction of eigenvectors, Figure 7A.PC1 collective motions extracted for the predominant eigenvectors of the using principal component in the studied complexes, Figure 7B.The scatter plots of the complexes in Figure, indicate the undergoing overall motions of the protein upon binding in terms of correlated and anti-correlated movements, Figure 7A and 7B.

Figure
Figure 3

3 A
Figure 3

Figure 7 A
Figure 7

Table 1 :
The best three compounds generated using XP docking with their corresponding docking scores and the interacting binding residues are presented.Catalytic dyad residues are shown in bold.

Table 2 :
In-silico ADMET screening of the selected compounds.
1 .Complex DE vdw DE elec DG gas DG polar DG nopolar DG solvation DG bind