In Silico Screening of Semi-Synthesized Compounds as Potential Inhibitors for SARS-CoV-2 Papain-like Protease: Pharmacophoric Features, Molecular Docking, ADMET, Toxicity and DFT Studies

Papain-like protease is an essential enzyme in the proteolytic processing required for the replication of SARS-CoV-2. Accordingly, such an enzyme is an important target for the development of anti-SARS-CoV-2 agents which may reduce the mortality associated with outbreaks of SARS-CoV-2. A set of 69 semi-synthesized molecules that exhibited the structural features of SARS-CoV-2 papain-like protease inhibitors (PLPI) were docked against the coronavirus papain-like protease (PLpro) enzyme (PDB ID: (4OW0). Docking studies showed that derivatives 34 and 58 were better than the co-crystallized ligand while derivatives 17, 28, 31, 40, 41, 43, 47, 54, and 65 exhibited good binding modes and binding free energies. The pharmacokinetic profiling study was conducted according to the four principles of the Lipinski rules and excluded derivative 31. Furthermore, ADMET and toxicity studies showed that derivatives 28, 34, and 47 have the potential to be drugs and have been demonstrated as safe when assessed via seven toxicity models. Finally, comparing the molecular orbital energies and the molecular electrostatic potential maps of 28, 34, and 47 against the co-crystallized ligand in a DFT study indicated that 28 is the most promising candidate to interact with the target receptor (PLpro).

These approaches have been practiced fruitfully and frequently in several scientific studies to find a treatment against COVID-19 with the advantage of taking less effort, time, and cost [17][18][19][20][21].

Docking Studies
MOE software was used to conduct docking studies (Supplementary Materials) on the investigated derivatives, with co-crystallized ligand S88 as a reference. The study aimed at getting a deeper insight into the binding modes of the examined semi-synthesized molecules in the active site of the coronavirus papain-like protease (PLpro) enzyme (PDB ID: (4OW0)). The docking method was validated through redocking of the co-crystallized ligand in the enzyme active site. The protocol's applicability was confirmed through the demonstration of small RMSD (0.54 Å) between the co-crystallized pose and the re-docked one ( Figure 4). In this study, we relied on the corrected mode of binding of the examined semisynthesized molecules and S88 as well as the values of the binding free energy (∆G) between them. Table 1 illustrates the calculated ∆G of the tested semi-synthesized molecules and (S88) against the coronavirus papain-like protease enzyme. The semi-synthesized molecules 34 and 58 showed affinity values of −8.97 and −8.65, respectively, that were higher than that of the redocked ligand S88 (−8.59 kcal/mol). Moreover, the semisynthesized molecules 17, 28, 31, 40, 41, 43, 47, 54 and 65 revealed binding energy scores ranging from −8.33 to −8.57 kcal/mol, which were highly close to the redocked ligand S88. On the other hand, the other semi-synthesized molecules demonstrated affinity values lower than S88. The proposed binding mode of the redocked ligand S88 revealed an affinity value of −8.59 kcal/mol. This high binding affinity is probably attributed to the formation of two hydrogen-bonding interactions. One was formed between the N-H group of the amide moiety and Tyr269 while the other was formed between the nitrogen atom of the pyridine ring and Asp165. Additionally, the naphthyl moiety formed four hydrophobic interactions with Tyr269 and Pro249. These results were found to be consistent with the reported data [45] ( Figure 5).
The docking simulation of compound 28 revealed that it has a good fitting into the enzyme active site with a docking score of −8.48 kcal/mol. The oxygen of the carbonyl group formed one hydrogen bond with the essential amino acid Tyr269. Additionally, the NH group of the pyrrole ring formed one hydrogen bond with Ala247. The 2,3,4,5tetrahydro-1H-pyrido [4,3-b]indole moiety formed Tyr265, Asp165, and Pro248 ( Figure 6).
As illustrated in Figure 7, compound 34 possessed a significant potential binding affinity (∆G = −8.97 kcal/mol) into the papain-like protease active site. This high binding affinity, which is higher than ligand S88, presumably attributed to the formation of one hydrogen bond interaction with Arg167. In addition, the 1H-indole formed two hydrophobic interactions with Lys158 and Leu163.
Investigation of the top docking poses of compounds 47 and 54 (affinity values of −8.57 and 8.33 kcal/mol) respectively, demonstrated that compound 47 formed one hydrogen bond interaction with Tyr269. In addition, it formed two hydrophobic interactions with Pro249 and Lys158. Compound 54 formed two hydrogen bonds with the essential amino acid Tyr269 and Tyr265. In addition, it formed two hydrophobic interactions with Lys158 and Arg167 in the active site of PLpro. Figures 8 and 9.
The binding mode of compound 58 (affinity value of −8.65 kcal/mol) was better than ligand S88. In detail, the amide moiety formed one hydrogen bond with fundamental amino acid Asp165, and NH of the pyrrole ring formed another hydrogen bond with Glu168. Furthermore, it formed aromatic stacking interactions (4 pi-cation bonds) with Tyr269, Tyr265, and Arg167 ( Figure 10).

Pharmacokinetic Profiling Study
An in silico computational evaluation of the physicochemical properties and profiling pharmacokinetics for the most active eleven semi-synthesized molecules, with ligand S88 as a reference compound, were conducted. The oral absorption of a drug is more likely to be better if the molecule fulfills at least three of the four principles of the Lipinski rules, listed below: (1) H bond donors (OH, NH, and SH) ≤ 5; (2) H bond acceptors (N, O, and S atoms) ≤ 10; (3) molecular weight < 500; (4) log P < 5. The bioavailability of compounds that violate more than one of these requirements is unlikely to be high. Moreover, reduced molecular flexibility, as measured by the number of rotatable bonds, and low polar surface area are found to be important predictors of good oral bioavailability [50]. Compounds with 10 or fewer rotatable bonds and a polar surface area of 140 Å or less have a high probability of good oral bioavailability [50,51]. The results, given in Table 2, revealed that all tested semi-synthesized molecules and reference ligand S88 showed no violation of Lipinski's rule except compound 31 (the Log P of compounds 31 was anticipated to be more than 5).

ADMET Studies
S88 and favipiravir were used as reference drugs in ADMET studies for the most active eleven semi-synthesized molecules using Discovery studio 4.0 software. ADMET studies include many descriptors. The predicted descriptors are listed in Table 3. All tested semi-synthesized molecules and favipiravir showed BBB penetration levels ranging from medium to low except compound 31, which displayed a very low BBB penetration level, and ligand S88 showed a high BBB penetration level. All semi-synthesized molecules, favipiravir, and ligand S88 have good absorption behavior except compound 31, which is expected to have a moderate absorption level. Moreover, the solubility level of the semi-synthesized molecules is projected to be better than or even comparable to that of the S88, which showed a low solubility level, except compound 31 that showed a very low solubility level. On the other hand, favipiravir demonstrated an optimal solubility level. All examined semi-synthesized molecules and favipiravir were predicted to be noninhibitors of CYP2D6 except compounds 31, 34, 47, and S88. Hepatotoxicity predictions found that all of the tested compounds and ligand S88 are predicted to be non-toxic except compounds 17, 31, 41, 43 and favipiravir, which have unfavorable hepatotoxic effects. All tested semi-synthesized molecules and S88 were expected to bind to plasma proteins more than 90% except compounds 28, 40, 43, 54, and favipiravir ( Figure 11).
As shown in Table 4, most of the examined semi-synthesized molecules have low toxicity. All the tested semi-synthesized molecules are non-carcinogens except 54, 58, and S88, which were predicted to be carcinogens. All tested semi-synthesized molecules showed TD 50 values ranging from 0.31 to 1.86 mg/kg body weight/day, which were lower than S88 (2.54 mg/kg body weight/day), except compounds 17, 31, 40, and 47 that showed TD 50 values of 5.57, 5.32, 10.31 and 14.54 mg/kg body weight/day, respectively, which were higher than S88. All the investigated semi-synthesized molecules revealed a maximum tolerated dose with a range of 0.045 to 0.122 g/kg body weight that was lower than S88 (0.124 g/kg body weight) except compounds 16, 28, 31, and 34, which demonstrated maximum tolerated doses of 0.142, 0.144, 0.369 and 0.328 g/kg body weight, respectively, which were higher than S88. All tested semi-synthesized molecules showed oral LD 50 values ranging from 2.97 to 32.39 mg/kg body weight/day that were higher than the LD 50 value of S88 (1.229 mg/kg body weight/day), except compound 31, which revealed an oral LD 50 value of 0.251 mg/kg body weight/day, which was lower than S88. Semi-synthesized molecules 54, 58, and 65 showed LOAEL values of 0.015, 0.001, and 0.018 g/kg body weight, respectively. These values were lower than S88 (0.035 g/kg body weight), while other semi-synthesized molecules revealed rat chronic LOAEL values ranging from 0.039 to 0.539 g/kg body weight, which were higher than S88. Moreover, all the tested semi-synthesized molecules and S88 were predicted to be mild irritants against the ocular irritancy model. On the other hand, the examined semi-synthesized molecules and S88 were expected to be non-irritant against the skin irritancy model except compound 58, which was anticipated to be a mild skin irritant. Figure 11. The expected ADMET study of the most potent semi-synthesized molecules.

Molecular Orbital Analysis
The total energies of 28, 34, 47, and S88 were −1422.912, −1285.184, −1252.334, and 1242.947kcal/mol, respectively. These results indicated that compound 28 has higher total energy than 34 and 47 and is expected to have a more efficient interaction with PLpro. Accordingly, compound 28 can bind more efficiently with PLpro than 34 and 47. Furthermore, the dipole moment values of 28, 34, 47, and S88 were 2.790, 1.558, 2.249, and 3.542, respectively (Table 5). These values indicate that 28 has a higher dipole moment than compounds 34 and 47. Based on these findings, it was expected that compound 28 can easily form hydrogen bonds and non-bonded interactions with PLpro, which, consequently, leads to an increased binding affinity with the target receptor during SARS-CoV-2 inhibition. Therefore, compound 28 is considered the most promising candidate to interact with the target receptor. As reported, HOMO and LUMO have a key role in chemical stability and reactivity [67]. Compound 28 had a gap energy value of 0.134 kcal/mol, which is higher than that of compounds 34 (0.099 kcal/mol) and 47 (0.097kcal/mol). The increased gap energy of compound 28 indicates the higher stability of this compound. Figure 12 showed the spatial distribution of molecular orbitals for the tested compounds.

Molecular Electrostatic Potential Maps (MEP)
MEP demonstrates the total electrostatic potential of a molecule in three dimensions depending on its partial charges, electronegativity, and chemical reactivity [68]. Identifying the electrostatic potential will help in the understanding of the drug's binding mode against a PLpro [69].
MEP displays the electronegative atoms (negative values) in red. Electronegative atoms act as hydrogen bonding acceptors. On the other hand, it displays electron-poor atoms (positive value) in blue. Electron-poor atoms act as hydrogen bonding donors. It displays the neutral atoms (zero values) in a green to yellow color. Neutral atoms can form πand other types of hydrophobic interactions. Such information facilitates the prediction of the chemical reaction and the binding mode with the biological target [70]. Compound 28 showed five red patches and two blue patches, which can form hydrogen bond acceptors and hydrogen bond donors, respectively. The aromatic moieties showed yellow patches, which can form hydrophobic interactions with hydrophobic amino acid residues (Figures 12 and 13).
Compounds 34 and 47 showed four red patches, which can form hydrogen bond acceptors. Compound 34 showed three red patches and two blue patches. The aromatic moieties of these compounds showed yellow patches which can form hydrophobic interactions with hydrophobic amino acid residues (Figures 12 and 13).
As compound 28 showed five red patches, this explains its high binding energy (−8.48 kcal/mol) and ability to form two hydrogen bonds in the docking procedure. The yellow patches on aromatic moieties facilitated the hydrophobic interaction with the target receptor. Compounds 34 and 47 showed four red patches in MEP, which clarified the formation of two and three hydrogen bonds, respectively. In addition, these compounds showed high binding energies of −8.97 and −8.57 kcal/mol, respectively.

Molecular Electrostatic Potential Maps (MEP)
MEP demonstrates the total electrostatic potential of a molecule in three dimensions depending on its partial charges, electronegativity, and chemical reactivity [68]. Identifying the electrostatic potential will help in the understanding of the drug's binding mode against a PLpro [69].
MEP displays the electronegative atoms (negative values) in red. Electronegative atoms act as hydrogen bonding acceptors. On the other hand, it displays electron-poor atoms (positive value) in blue. Electron-poor atoms act as hydrogen bonding donors. It displays the neutral atoms (zero values) in a green to yellow color. Neutral atoms can form π-and other types of hydrophobic interactions. Such information facilitates the prediction of the chemical reaction and the binding mode with the biological target [70].
Compound 28 showed five red patches and two blue patches, which can form hydrogen bond acceptors and hydrogen bond donors, respectively. The aromatic moieties showed yellow patches, which can form hydrophobic interactions with hydrophobic As compound 28 showed five red patches, this explains its high binding energy (−8.48 kcal/mol) and ability to form two hydrogen bonds in the docking procedure. The yellow patches on aromatic moieties facilitated the hydrophobic interaction with the target receptor. Compounds 34 and 47 showed four red patches in MEP, which clarified the formation of two and three hydrogen bonds, respectively. In addition, these compounds showed high binding energies of −8.97 and −8.57 kcal/mol, respectively.

Conclusions
A set of 69 semi-synthesized molecules that exhibited the structural features of PLpro inhibitors (PLPI) were screened in silico to select the most potent inhibitor of PLpro enzyme (PDB ID: (4OW0). Docking studies showed that 11 molecules exhibited good binding modes and binding free energies. The pharmacokinetic profiling study excluded an unsuitable compound. Furthermore, ADMET and toxicity studies favored three molecules. Finally, a DFT study has been carried out and indicated that N-(3,4-dimethoxyphenethyl)-4-oxo-4-(1,3,4,5-tetrahydro-2H-pyrido [4,3-b]indol-2-yl)butanamide (28) is the most promising PLpro inhibitor. Further work has to be conducted to build on the presented results in the hopes of finding a cure.

Conclusions
A set of 69 semi-synthesized molecules that exhibited the structural features of PLpro inhibitors (PLPI) were screened in silico to select the most potent inhibitor of PLpro enzyme (PDB ID: (4OW0). Docking studies showed that 11 molecules exhibited good binding modes and binding free energies. The pharmacokinetic profiling study excluded an unsuitable compound. Furthermore, ADMET and toxicity studies favored three molecules. Finally, a DFT study has been carried out and indicated that N-(3,4-dimethoxyphenethyl)-4-oxo-4-(1,3,4,5-tetrahydro-2H-pyrido [4,3-b]indol-2-yl)butanamide (28) is the most promising PLpro inhibitor. Further work has to be conducted to build on the presented results in the hopes of finding a cure.

Docking Studies
Crystal structure of human coronavirus papain-like protease inhibitor [PDB ID: 4OW0, resolution: 2.10 Å] was obtained from Protein Data Bank. The docking investigation was accomplished using MOE2014 software. At first, the crystal structure of SARS-CoV-2 helicase was prepared by removing water molecules. Only one chain was retained beside the co-crystallized ligand (S88). Then, the selected chain was protonated and subjected to minimization of energy process. Next, the active site of the target protein was defined.
Structures of the tested compounds and the co-crystallized ligand were drawn using ChemBioDraw Ultra 14.0 and saved as MDL-SD format. Such file was opened using MOE to display the 3D structures which were protonated and subjected to energy minimization. Formerly, validation of the docking process was performed by docking the co-crystallized ligand against the isolated pocket of active site. The produced RMSD value indicated the validity of process. Finally, docking of the tested compounds was done through the dock option inserted in computer window. For each docked molecule, 30 docked poses were produced using ASE for scoring function and force field for refinement. The retain was kept at 30. The crystal parameters were adjusted at default values (Coordinates: Normal, Lattice Style: the same, Lattice: (1)). The results of the docking process were then visualized using Discovery Studio 4.0 software [71].

Pharmacokinetic Profiling
Pharmacokinetic profile of the compounds was determined using Discovery studio 4.0. [72].

ADMET Analysis
ADMET descriptors (absorption, distribution, metabolism, excretion and toxicity) of the compounds were determined using Discovery studio 4.0. At first, the CHARMM force field was applied, then the tested compounds were prepared and minimized according to the preparation of small molecule protocol. Then ADMET descriptors protocol was applied to carry out these studies [73].

Toxicity Studies
The toxicity parameters of the tested compounds were calculated using Discovery studio 4.0. Indinavir was used as a reference drug. At first, the CHARMM force field was applied then the compounds were prepared and minimized according to the preparation of small molecule protocol. Then different parameters were calculated from the toxicity prediction (extensible) protocol [73][74][75].

DFT Studies
The DFT parameters (total energy, binding energy, HOMO, LUMO, gap energy, dipole moment, and electrostatic potential) were calculated using Discovery studio software. The tested compounds were prepared using prepare ligand protocol. Then, the prepared compounds were subjected to DFT calculation protocol using the default option [76].
Supplementary Materials: The following are available online, Table S1: Detailed toxicity report, in addition to the method (Docking studies, ADMET studies, Toxicity studies and DFT studies).