Molecular Docking and Dynamics Investigations for Identifying Potential Inhibitors of the 3-Chymotrypsin-like Protease of SARS-CoV-2: Repurposing of Approved Pyrimidonic Pharmaceuticals for COVID-19 Treatment

This study demonstrates the inhibitory effect of 42 pyrimidonic pharmaceuticals (PPs) on the 3-chymotrypsin-like protease of SARS-CoV-2 (3CLpro) through molecular docking, molecular dynamics simulations, and free binding energies by means of molecular mechanics–Poisson Boltzmann surface area (MM-PBSA) and molecular mechanics–generalized Born surface area (MM-GBSA). Of these tested PPs, 11 drugs approved by the US Food and Drug Administration showed an excellent binding affinity to the catalytic residues of 3CLpro of His41 and Cys145: uracil mustard, cytarabine, floxuridine, trifluridine, stavudine, lamivudine, zalcitabine, telbivudine, tipiracil, citicoline, and uridine triacetate. Their percentage of residues involved in binding at the active sites ranged from 56 to 100, and their binding affinities were in the range from −4.6 ± 0.14 to −7.0 ± 0.19 kcal/mol. The molecular dynamics as determined by a 200 ns simulation run of solvated docked complexes confirmed the stability of PP conformations that bound to the catalytic dyad and the active sites of 3CLpro. The free energy of binding also demonstrates the stability of the PP–3CLpro complexes. Citicoline and uridine triacetate showed free binding energies of −25.53 and −7.07 kcal/mol, respectively. Therefore, I recommend that they be repurposed for the fight against COVID-19, following proper experimental and clinical validation.


Introduction
Over a year has passed since the COVID-19 pandemic began. Some vaccines, such as those by Pfizer and Moderna, and some drugs, such as remdesivir, have been approved for use in therapy. The efforts by governments, health organizations, and other sectors to stem the alarmingly increasing numbers of deaths and cases were unprecedented [1][2][3][4][5][6]. However, SARS-CoV-2 continues to threaten the world, with over four million deaths and 227 million cases as of 16 September 2021 (https://www.worldometers.info/coronavirus/ accessed on 27 November 2021). COVID-19 was declared a pandemic by the World Health Organization on 11 March 2020. Today, the new SARS-CoV-2 virus, the causative agent of COVID-19, has been detected in almost every country on the planet [5,[7][8][9].
Coronaviruses are positive-stranded RNA viruses with the largest viral genomes ever known, ranging from 16 to 32 kb. The 3-chymotrypsin-like protease (3CL pro ) produced by SARS-CoV-2 is a cysteine protease encoded as nonstructural protein 3 in the polyprotein. 3CL pro is responsible for the cleavage of 11 specific sites of polyproteins (pp1a, pp1ab) produced by the 229E gene. These polyproteins are involved in the production of a functional polypeptide essential for viral replication and transcription. Further, the specificity of 3CL pro is dissimilar to that of human host-cell protease. Thus, 3CL pro has become the focus of drug repurposing and development programs to combat the COVID-19 pandemic [10][11][12][13]. ther, the specificity of 3CL pro is dissimilar to that of human host-cell protease. Thus, 3CL pro has become the focus of drug repurposing and development programs to combat the COVID-19 pandemic [10][11][12][13]. Recent and ongoing research has reported that some pharmaceutical, synthetic, and natural products can act as 3CL pro inhibitors or against SARS-CoV-2 in general. These include selenium-containing heterocyclic compounds, chloroquine phosphate, indinavir, darunavir, lopinavir, eravacycline, naproxen, salix cortex, antioxidants, chiral phytochemicals from Opuntia ficus-indica, elbasvir, valrubicin, favipiravir isoflavone, and myricitrin [6,[14][15][16][17][18][19][20][21][22]. Although some of these have entered human clinical trials or were even approved, more studies are still needed. The importance of the pyridone ring was highlighted in synthetic materials and drugs containing pyridone [11,23]. The pyrimidone ring has the exact shape of pyridine but is more functionalized and electrondeficient. Herein, we screened the inhibitory activity of 42 approved pyrimidonic pharmaceuticals (PPs) against 3CL pro using a combination of molecular docking analyses, molecular dynamics simulations, and calculations of the MM-PBSA and MM-GBSA binding free energies. The sites of action of active inhibitors were investigated, discussed, and explored.

The Pyrimidonic Pharmaceuticals (PPs)
The PPs were selected using the search engine of the drug bank database. The search uncovered 46 PPs; the pharmaceuticals containing caffeine were entirely excluded as all except enprofylline have previously been studied. In addition, the macropolymeric drug mipomersen, drugs composed of a mixture of drugs, and withdrawn drugs were not included in this study. The chosen drugs were classified into four categories according to their structures. 1PPs have only one heterocycle, 2aPPs have two, 2bPPs have two heterocycles with a pyrimidone ring having an extra carbonyl group, and 3PPs have three or more heterocycles (Table 1). ther, the specificity of 3CL pro is dissimilar to that of human host-cell protease. Thus, 3CL pro has become the focus of drug repurposing and development programs to combat the COVID-19 pandemic [10][11][12][13]. Recent and ongoing research has reported that some pharmaceutical, synthetic, and natural products can act as 3CL pro inhibitors or against SARS-CoV-2 in general. These include selenium-containing heterocyclic compounds, chloroquine phosphate, indinavir, darunavir, lopinavir, eravacycline, naproxen, salix cortex, antioxidants, chiral phytochemicals from Opuntia ficus-indica, elbasvir, valrubicin, favipiravir isoflavone, and myricitrin [6,[14][15][16][17][18][19][20][21][22]. Although some of these have entered human clinical trials or were even approved, more studies are still needed. The importance of the pyridone ring was highlighted in synthetic materials and drugs containing pyridone [11,23]. The pyrimidone ring has the exact shape of pyridine but is more functionalized and electrondeficient. Herein, we screened the inhibitory activity of 42 approved pyrimidonic pharmaceuticals (PPs) against 3CL pro using a combination of molecular docking analyses, molecular dynamics simulations, and calculations of the MM-PBSA and MM-GBSA binding free energies. The sites of action of active inhibitors were investigated, discussed, and explored.

The Pyrimidonic Pharmaceuticals (PPs)
The PPs were selected using the search engine of the drug bank database. The search uncovered 46 PPs; the pharmaceuticals containing caffeine were entirely excluded as all except enprofylline have previously been studied. In addition, the macropolymeric drug mipomersen, drugs composed of a mixture of drugs, and withdrawn drugs were not included in this study. The chosen drugs were classified into four categories according to their structures. 1PPs have only one heterocycle, 2aPPs have two, 2bPPs have two heterocycles with a pyrimidone ring having an extra carbonyl group, and 3PPs have three or more heterocycles (Table 1). ther, the specificity of 3CL pro is dissimilar to that of human host-cell protease. Thus, 3CL pro has become the focus of drug repurposing and development programs to combat the COVID-19 pandemic [10][11][12][13]. Recent and ongoing research has reported that some pharmaceutical, synthetic, and natural products can act as 3CL pro inhibitors or against SARS-CoV-2 in general. These include selenium-containing heterocyclic compounds, chloroquine phosphate, indinavir, darunavir, lopinavir, eravacycline, naproxen, salix cortex, antioxidants, chiral phytochemicals from Opuntia ficus-indica, elbasvir, valrubicin, favipiravir isoflavone, and myricitrin [6,[14][15][16][17][18][19][20][21][22]. Although some of these have entered human clinical trials or were even approved, more studies are still needed. The importance of the pyridone ring was highlighted in synthetic materials and drugs containing pyridone [11,23]. The pyrimidone ring has the exact shape of pyridine but is more functionalized and electrondeficient. Herein, we screened the inhibitory activity of 42 approved pyrimidonic pharmaceuticals (PPs) against 3CL pro using a combination of molecular docking analyses, molecular dynamics simulations, and calculations of the MM-PBSA and MM-GBSA binding free energies. The sites of action of active inhibitors were investigated, discussed, and explored.

The Pyrimidonic Pharmaceuticals (PPs)
The PPs were selected using the search engine of the drug bank database. The search uncovered 46 PPs; the pharmaceuticals containing caffeine were entirely excluded as all except enprofylline have previously been studied. In addition, the macropolymeric drug mipomersen, drugs composed of a mixture of drugs, and withdrawn drugs were not included in this study. The chosen drugs were classified into four categories according to their structures. 1PPs have only one heterocycle, 2aPPs have two, 2bPPs have two heterocycles with a pyrimidone ring having an extra carbonyl group, and 3PPs have three or more heterocycles (Table 1). ther, the specificity of 3CL pro is dissimilar to that of human host-cell protease. Thus, 3CL pro has become the focus of drug repurposing and development programs to combat the COVID-19 pandemic [10][11][12][13]. Recent and ongoing research has reported that some pharmaceutical, synthetic, and natural products can act as 3CL pro inhibitors or against SARS-CoV-2 in general. These include selenium-containing heterocyclic compounds, chloroquine phosphate, indinavir, darunavir, lopinavir, eravacycline, naproxen, salix cortex, antioxidants, chiral phytochemicals from Opuntia ficus-indica, elbasvir, valrubicin, favipiravir isoflavone, and myricitrin [6,[14][15][16][17][18][19][20][21][22]. Although some of these have entered human clinical trials or were even approved, more studies are still needed. The importance of the pyridone ring was highlighted in synthetic materials and drugs containing pyridone [11,23]. The pyrimidone ring has the exact shape of pyridine but is more functionalized and electrondeficient. Herein, we screened the inhibitory activity of 42 approved pyrimidonic pharmaceuticals (PPs) against 3CL pro using a combination of molecular docking analyses, molecular dynamics simulations, and calculations of the MM-PBSA and MM-GBSA binding free energies. The sites of action of active inhibitors were investigated, discussed, and explored.

The Pyrimidonic Pharmaceuticals (PPs)
The PPs were selected using the search engine of the drug bank database. The search uncovered 46 PPs; the pharmaceuticals containing caffeine were entirely excluded as all except enprofylline have previously been studied. In addition, the macropolymeric drug mipomersen, drugs composed of a mixture of drugs, and withdrawn drugs were not included in this study. The chosen drugs were classified into four categories according to their structures. 1PPs have only one heterocycle, 2aPPs have two, 2bPPs have two heterocycles with a pyrimidone ring having an extra carbonyl group, and 3PPs have three or more heterocycles (Table 1). ther, the specificity of 3CL pro is dissimilar to that of human host-cell protease. Thus, 3CL pro has become the focus of drug repurposing and development programs to combat the COVID-19 pandemic [10][11][12][13]. Recent and ongoing research has reported that some pharmaceutical, synthetic, and natural products can act as 3CL pro inhibitors or against SARS-CoV-2 in general. These include selenium-containing heterocyclic compounds, chloroquine phosphate, indinavir, darunavir, lopinavir, eravacycline, naproxen, salix cortex, antioxidants, chiral phytochemicals from Opuntia ficus-indica, elbasvir, valrubicin, favipiravir isoflavone, and myricitrin [6,[14][15][16][17][18][19][20][21][22]. Although some of these have entered human clinical trials or were even approved, more studies are still needed. The importance of the pyridone ring was highlighted in synthetic materials and drugs containing pyridone [11,23]. The pyrimidone ring has the exact shape of pyridine but is more functionalized and electrondeficient. Herein, we screened the inhibitory activity of 42 approved pyrimidonic pharmaceuticals (PPs) against 3CL pro using a combination of molecular docking analyses, molecular dynamics simulations, and calculations of the MM-PBSA and MM-GBSA binding free energies. The sites of action of active inhibitors were investigated, discussed, and explored.

The Pyrimidonic Pharmaceuticals (PPs)
The PPs were selected using the search engine of the drug bank database. The search uncovered 46 PPs; the pharmaceuticals containing caffeine were entirely excluded as all except enprofylline have previously been studied. In addition, the macropolymeric drug mipomersen, drugs composed of a mixture of drugs, and withdrawn drugs were not included in this study. The chosen drugs were classified into four categories according to their structures. 1PPs have only one heterocycle, 2aPPs have two, 2bPPs have two heterocycles with a pyrimidone ring having an extra carbonyl group, and 3PPs have three or more heterocycles (Table 1). ther, the specificity of 3CL pro is dissimilar to that of human host-cell protease. Thus, 3CL pro has become the focus of drug repurposing and development programs to combat the COVID-19 pandemic [10][11][12][13]. Recent and ongoing research has reported that some pharmaceutical, synthetic, and natural products can act as 3CL pro inhibitors or against SARS-CoV-2 in general. These include selenium-containing heterocyclic compounds, chloroquine phosphate, indinavir, darunavir, lopinavir, eravacycline, naproxen, salix cortex, antioxidants, chiral phytochemicals from Opuntia ficus-indica, elbasvir, valrubicin, favipiravir isoflavone, and myricitrin [6,[14][15][16][17][18][19][20][21][22]. Although some of these have entered human clinical trials or were even approved, more studies are still needed. The importance of the pyridone ring was highlighted in synthetic materials and drugs containing pyridone [11,23]. The pyrimidone ring has the exact shape of pyridine but is more functionalized and electrondeficient. Herein, we screened the inhibitory activity of 42 approved pyrimidonic pharmaceuticals (PPs) against 3CL pro using a combination of molecular docking analyses, molecular dynamics simulations, and calculations of the MM-PBSA and MM-GBSA binding free energies. The sites of action of active inhibitors were investigated, discussed, and explored.

The Pyrimidonic Pharmaceuticals (PPs)
The PPs were selected using the search engine of the drug bank database. The search uncovered 46 PPs; the pharmaceuticals containing caffeine were entirely excluded as all except enprofylline have previously been studied. In addition, the macropolymeric drug mipomersen, drugs composed of a mixture of drugs, and withdrawn drugs were not included in this study. The chosen drugs were classified into four categories according to their structures. 1PPs have only one heterocycle, 2aPPs have two, 2bPPs have two heterocycles with a pyrimidone ring having an extra carbonyl group, and 3PPs have three or more heterocycles (Table 1). ther, the specificity of 3CL pro is dissimilar to that of human host-cell protease. Thus, 3CL pro has become the focus of drug repurposing and development programs to combat the COVID-19 pandemic [10][11][12][13]. Recent and ongoing research has reported that some pharmaceutical, synthetic, and natural products can act as 3CL pro inhibitors or against SARS-CoV-2 in general. These include selenium-containing heterocyclic compounds, chloroquine phosphate, indinavir, darunavir, lopinavir, eravacycline, naproxen, salix cortex, antioxidants, chiral phytochemicals from Opuntia ficus-indica, elbasvir, valrubicin, favipiravir isoflavone, and myricitrin [6,[14][15][16][17][18][19][20][21][22]. Although some of these have entered human clinical trials or were even approved, more studies are still needed. The importance of the pyridone ring was highlighted in synthetic materials and drugs containing pyridone [11,23]. The pyrimidone ring has the exact shape of pyridine but is more functionalized and electrondeficient. Herein, we screened the inhibitory activity of 42 approved pyrimidonic pharmaceuticals (PPs) against 3CL pro using a combination of molecular docking analyses, molecular dynamics simulations, and calculations of the MM-PBSA and MM-GBSA binding free energies. The sites of action of active inhibitors were investigated, discussed, and explored.

The Pyrimidonic Pharmaceuticals (PPs)
The PPs were selected using the search engine of the drug bank database. The search uncovered 46 PPs; the pharmaceuticals containing caffeine were entirely excluded as all except enprofylline have previously been studied. In addition, the macropolymeric drug mipomersen, drugs composed of a mixture of drugs, and withdrawn drugs were not included in this study. The chosen drugs were classified into four categories according to their structures. 1PPs have only one heterocycle, 2aPPs have two, 2bPPs have two heterocycles with a pyrimidone ring having an extra carbonyl group, and 3PPs have three or more heterocycles (Table 1). ther, the specificity of 3CL pro is dissimilar to that of human host-cell protease. Thus, 3CL pro has become the focus of drug repurposing and development programs to combat the COVID-19 pandemic [10][11][12][13]. Recent and ongoing research has reported that some pharmaceutical, synthetic, and natural products can act as 3CL pro inhibitors or against SARS-CoV-2 in general. These include selenium-containing heterocyclic compounds, chloroquine phosphate, indinavir, darunavir, lopinavir, eravacycline, naproxen, salix cortex, antioxidants, chiral phytochemicals from Opuntia ficus-indica, elbasvir, valrubicin, favipiravir isoflavone, and myricitrin [6,[14][15][16][17][18][19][20][21][22]. Although some of these have entered human clinical trials or were even approved, more studies are still needed. The importance of the pyridone ring was highlighted in synthetic materials and drugs containing pyridone [11,23]. The pyrimidone ring has the exact shape of pyridine but is more functionalized and electrondeficient. Herein, we screened the inhibitory activity of 42 approved pyrimidonic pharmaceuticals (PPs) against 3CL pro using a combination of molecular docking analyses, molecular dynamics simulations, and calculations of the MM-PBSA and MM-GBSA binding free energies. The sites of action of active inhibitors were investigated, discussed, and explored.

The Pyrimidonic Pharmaceuticals (PPs)
The PPs were selected using the search engine of the drug bank database. The search uncovered 46 PPs; the pharmaceuticals containing caffeine were entirely excluded as all except enprofylline have previously been studied. In addition, the macropolymeric drug mipomersen, drugs composed of a mixture of drugs, and withdrawn drugs were not included in this study. The chosen drugs were classified into four categories according to their structures. 1PPs have only one heterocycle, 2aPPs have two, 2bPPs have two heterocycles with a pyrimidone ring having an extra carbonyl group, and 3PPs have three or more heterocycles (Table 1).

Generation and Energy Minimization of the PPs and 3CL pro
The 3D structures of the selected PPs were downloaded from the PubChem website as SDF files; their energy was minimized for 10,000 steepest descent steps at 5000 conjugate gradient steps using antechamber plugin UCSF Chimera [24,25]. For alogliptin, the 3D structure was obtained by utilizing OpenBabel converter tools and ChemSkech [26]. The crystal structure of SARS-CoV-2 3CL pro was obtained from the Protein Data Bank database website (PDB ID: 6Y2E). For analysis, water was removed from the 3CL pro structure, and the energy was then minimized for 1000 steepest descent steps at 20 conjugate gradient steps.

Molecular Docking
Blind molecular docking experiments were performed using the AutoDock Vina tool implemented with the interactive visualization and analysis program UCSF Chimera. The default parameter values were adopted with a grid box (−15 × −25 × 15) Å, centered at (35,65, 65) Å. The predicted affinity values of the score were observed using the View Dock tool. The binding between ligands and 3CL pro active sites and the images were processed and visualized using UCSF Chimera [24][25][26][27][28].

Molecular Dynamics Simulations
MD simulations were performed as previously described [29]. The PP ligands were separated from the docked complexes using UCSF Chimera. The missed hydrogens were added and saved as PDB files using AMBER's large-structure serial numbering. Topology files and parameters of the receptor and the ligands were made using leap and antechamber of Amber Tools 21 [30,30], utilizing Amber force fields of GAFF2 [31] and ff14SB [32] to assign inhibitors and 3CL pro structure, respectively. The systems were solvated with TIP3P water molecules [33] and were neutralized via sodium ions. Subsequently, molecular dynamics (MD) simulations were performed by means of the Nanoscale Molecular Dynamics (NAMD) Simulation 2.6 program [34]. Each system was

Generation and Energy Minimization of the PPs and 3CL pro
The 3D structures of the selected PPs were downloaded from the PubChem website as SDF files; their energy was minimized for 10,000 steepest descent steps at 5000 conjugate gradient steps using antechamber plugin UCSF Chimera [24,25]. For alogliptin, the 3D structure was obtained by utilizing OpenBabel converter tools and ChemSkech [26]. The crystal structure of SARS-CoV-2 3CL pro was obtained from the Protein Data Bank database website (PDB ID: 6Y2E). For analysis, water was removed from the 3CL pro structure, and the energy was then minimized for 1000 steepest descent steps at 20 conjugate gradient steps.

Molecular Docking
Blind molecular docking experiments were performed using the AutoDock Vina tool implemented with the interactive visualization and analysis program UCSF Chimera. The default parameter values were adopted with a grid box (−15 × −25 × 15) Å, centered at (35,65, 65) Å. The predicted affinity values of the score were observed using the View Dock tool. The binding between ligands and 3CL pro active sites and the images were processed and visualized using UCSF Chimera [24][25][26][27][28].

Molecular Dynamics Simulations
MD simulations were performed as previously described [29]. The PP ligands were separated from the docked complexes using UCSF Chimera. The missed hydrogens were added and saved as PDB files using AMBER's large-structure serial numbering. Topology files and parameters of the receptor and the ligands were made using leap and antechamber of Amber Tools 21 [30,30], utilizing Amber force fields of GAFF2 [31] and ff14SB [32] to assign inhibitors and 3CL pro structure, respectively. The systems were solvated with TIP3P water molecules [33] and were neutralized via sodium ions. Subsequently, molecular dynamics (MD) simulations were performed by means of the Nanoscale Molecular Dynamics (NAMD) Simulation 2.6 program [34]. Each system was

Generation and Energy Minimization of the PPs and 3CL pro
The 3D structures of the selected PPs were downloaded from the PubChem website as SDF files; their energy was minimized for 10,000 steepest descent steps at 5000 conjugate gradient steps using antechamber plugin UCSF Chimera [24,25]. For alogliptin, the 3D structure was obtained by utilizing OpenBabel converter tools and ChemSkech [26]. The crystal structure of SARS-CoV-2 3CL pro was obtained from the Protein Data Bank database website (PDB ID: 6Y2E). For analysis, water was removed from the 3CL pro structure, and the energy was then minimized for 1000 steepest descent steps at 20 conjugate gradient steps.

Molecular Docking
Blind molecular docking experiments were performed using the AutoDock Vina tool implemented with the interactive visualization and analysis program UCSF Chimera. The default parameter values were adopted with a grid box (−15 × −25 × 15) Å, centered at (35,65, 65) Å. The predicted affinity values of the score were observed using the View Dock tool. The binding between ligands and 3CL pro active sites and the images were processed and visualized using UCSF Chimera [24][25][26][27][28].

Molecular Dynamics Simulations
MD simulations were performed as previously described [29]. The PP ligands were separated from the docked complexes using UCSF Chimera. The missed hydrogens were added and saved as PDB files using AMBER's large-structure serial numbering. Topology files and parameters of the receptor and the ligands were made using leap and antechamber of Amber Tools 21 [30,30], utilizing Amber force fields of GAFF2 [31] and ff14SB [32] to assign inhibitors and 3CL pro structure, respectively. The systems were solvated with TIP3P water molecules [33] and were neutralized via sodium ions. Subsequently, molecular dynamics (MD) simulations were performed by means of the Nanoscale Molecular Dynamics (NAMD) Simulation 2.6 program [34]. Each system was

Generation and Energy Minimization of the PPs and 3CL pro
The 3D structures of the selected PPs were downloaded from the PubChem website as SDF files; their energy was minimized for 10,000 steepest descent steps at 5000 conjugate gradient steps using antechamber plugin UCSF Chimera [24,25]. For alogliptin, the 3D structure was obtained by utilizing OpenBabel converter tools and ChemSkech [26]. The crystal structure of SARS-CoV-2 3CL pro was obtained from the Protein Data Bank database website (PDB ID: 6Y2E). For analysis, water was removed from the 3CL pro structure, and the energy was then minimized for 1000 steepest descent steps at 20 conjugate gradient steps.

Molecular Docking
Blind molecular docking experiments were performed using the AutoDock Vina tool implemented with the interactive visualization and analysis program UCSF Chimera. The default parameter values were adopted with a grid box (−15 × −25 × 15) Å, centered at (35, 65, 65) Å. The predicted affinity values of the score were observed using the View Dock tool. The binding between ligands and 3CL pro active sites and the images were processed and visualized using UCSF Chimera [24][25][26][27][28].

Molecular Dynamics Simulations
MD simulations were performed as previously described [29]. The PP ligands were separated from the docked complexes using UCSF Chimera. The missed hydrogens were added and saved as PDB files using AMBER's large-structure serial numbering. Topology files and parameters of the receptor and the ligands were made using leap and antechamber of Amber Tools 21 [30,30], utilizing Amber force fields of GAFF2 [31] and ff14SB [32] to assign inhibitors and 3CL pro structure, respectively. The systems were solvated with TIP3P water molecules [33] and were neutralized via sodium ions. Subsequently, molecular dynamics (MD) simulations were performed by means of the Nanoscale Molecular Dynamics (NAMD) Simulation 2.6 program [34]. Each system was -

Generation and Energy Minimization of the PPs and 3CL pro
The 3D structures of the selected PPs were downloaded from the PubChem website as SDF files; their energy was minimized for 10,000 steepest descent steps at 5000 conjugate gradient steps using antechamber plugin UCSF Chimera [24,25]. For alogliptin, the 3D structure was obtained by utilizing OpenBabel converter tools and ChemSkech [26]. The crystal structure of SARS-CoV-2 3CL pro was obtained from the Protein Data Bank database website (PDB ID: 6Y2E). For analysis, water was removed from the 3CL pro structure, and the energy was then minimized for 1000 steepest descent steps at 20 conjugate gradient steps.

Molecular Docking
Blind molecular docking experiments were performed using the AutoDock Vina tool implemented with the interactive visualization and analysis program UCSF Chimera. The default parameter values were adopted with a grid box (−15 × −25 × 15) Å, centered at (35, 65, 65) Å. The predicted affinity values of the score were observed using the View Dock tool. The binding between ligands and 3CL pro active sites and the images were processed and visualized using UCSF Chimera [24][25][26][27][28].

Molecular Dynamics Simulations
MD simulations were performed as previously described [29]. The PP ligands were separated from the docked complexes using UCSF Chimera. The missed hydrogens were added and saved as PDB files using AMBER's large-structure serial numbering. Topology files and parameters of the receptor and the ligands were made using leap and antechamber of Amber Tools 21 [30,30], utilizing Amber force fields of GAFF2 [31] and ff14SB [32] to assign inhibitors and 3CL pro structure, respectively. The systems were solvated with TIP3P water molecules [33] and were neutralized via sodium ions. Subsequently, molecular dynamics (MD) simulations were performed by means of the Nanoscale Molecular Dynamics (NAMD) Simulation 2.6 program [34]. Each system was minimized for 1 ps at 273.15 K using the NVE ensemble. The temperature was gradually increased to 310 K using the NVT ensemble in a protocol consisting of 1600 minimization steps. Then, each system was minimized for 10 ps at 310 K followed by 200 ns of MD simulation control using the NVT ensemble at 310 K and a time step of 2 fs. In order to calculate electrostatic interactions, the particle mesh Ewald process and periodic boundary conditions were applied [35,36]. The root-mean-square fluctuation (RMSF) and the root-mean-square deviation (RMSD) for each system were obtained by analyzing the trajectory using the VMD 1.8 program [37].

The Binding Free Energies
The binding free energies of the PP-3CL pro complexes were calculated by means of molecular mechanics-Poisson Boltzmann surface area (MM-PBSA) and molecular mechanics-generalized Born surface area (MM-GBSA) using the MMPBSA.py module of Amber Tools 21 [38]. The MD simulation over 200 ns provided several conformations sampled after equilibrium, using the last frames to lessen the computational cost. CPPTRAJ was used to obtain the snapshots [39]. The conformational changes were evaluated through quasi-harmonic entropy approximation [40]. The free energy of the binding interaction between inhibitors and 3CL pro complexes can be obtained via the following equations: ∆G gas = E vdw + E elec (3) where ∆H represents enthalpy change, T∆S represents the entropic contribution, E vdw represents the van der Waals interaction energy, E ele represents the electrostatic interaction energy, ∆G sol represents the polar solvation energy, and E np represents the nonpolar solvation energy.

Results and Discussion
The results of the molecular docking are tabulated in Tables S1-S4. Figure 1 shows the catalytic dyad and the active sites of 3CL pro . The crucial residues HIS 41, GLY 143, SER 144, and CYS 145 forming the S1 site are shown in black. Then, PHE 140, LUE 141, ASN 142, HIS 163, GLU 166 (magenta), and the N-terminal amino acid residues (blue) are involved in the formation of the S1 subsite of the substrate-binding pocket. The MET 49, TYR 54, HIS 164, ASP 187, and ARG R188 residues form the S2 site (green). MET 165, LEU 167, GLN 189, THR 190, and GLN 192 comprise the S4 site (cyan). The SER 284, ALA 285, and LEU 286 residues (yellow) are a result of genetic mutation leading to an increase in the SASR-CoV-2 3CL pro activity of 3.6 fold over that of the 3CL pro predecessor of SARS-CoV [12,41]. the SASR-CoV-2 3CL pro activity of 3.6 fold over that of the 3CL pro predecessor of SARS-CoV [12,42].

Molecular Docking
The docked complexes of the top 11 candidates are depicted in Figure 2. Their binding affinities to the active sites of 3CL pro are shown in Table 2. The 3PPs showed significant interactions with the residues LEU 286, SER 284, and ALA 285, and a relatively lower interaction ratio to the catalytic dyad, in contrast to the other groups. Of the 3PPs, alogliptin and flavin mononucleotide were found to have the highest binding percentage with the catalytic dyad and to form hydrogen bonds with the S1 and S'1 sites. These were followed by riboflavin and sofosbuvir with an advantage in binding to the LEU 286 residue (Table S1). Flavin adenine dinucleotide showed excellent binding affinity to LEU 286 but not with the catalytic dyad. Zidovudine and gemcitabine demonstrated similar activity to alogliptin and flavin mononucleotide (Tables S2 and S3).
Among the 2bPPs, anti-hepatitis B infection telbivudine, anti-orotic aciduria uridine triacetate, and anticancer tipiracil were found to have the highest binding to 3CL pro active sites, followed by antimetabolite floxuridine, anti-herpesvirus trifluridine, and anti-HIV stavudine. Here, it is worth noting the importance of the molecular structure, as this set differed from the previous one by its increased ability to bind to the 3CL pro catalytic dyad. The 2aPPs showed similar activity to that of the 2bPPs. Anticancer cytarabine, anti- Color indicates the residues involved in the formation of the S1 site (shown in magenta), S1 site from the other promotor (blue), S2 site (green), S4 site (cyan), and S1' site (black), in addition to SER 284, ALA 285, and LEU 286 (yellow). (b) Only the catalytic dyad and GLU 166 residues.

Molecular Docking
The docked complexes of the top 11 candidates are depicted in Figure 2. Their binding affinities to the active sites of 3CL pro are shown in Table 2. The 3PPs showed significant interactions with the residues LEU 286, SER 284, and ALA 285, and a relatively lower interaction ratio to the catalytic dyad, in contrast to the other groups. Of the 3PPs, alogliptin and flavin mononucleotide were found to have the highest binding percentage with the catalytic dyad and to form hydrogen bonds with the S1 and S'1 sites. These were followed by riboflavin and sofosbuvir with an advantage in binding to the LEU 286 residue (Table S1). Flavin adenine dinucleotide showed excellent binding affinity to LEU 286 but not with the catalytic dyad. Zidovudine and gemcitabine demonstrated similar activity to alogliptin and flavin mononucleotide (Tables S2 and S3).
Among the 2bPPs, anti-hepatitis B infection telbivudine, anti-orotic aciduria uridine triacetate, and anticancer tipiracil were found to have the highest binding to 3CL pro active sites, followed by antimetabolite floxuridine, anti-herpesvirus trifluridine, and anti-HIV stavudine. Here, it is worth noting the importance of the molecular structure, as this set differed from the previous one by its increased ability to bind to the 3CL pro catalytic dyad. The 2aPPs showed similar activity to that of the 2bPPs. Anticancer cytarabine, antiglaucoma citicoline, and anti-HIV drugs lamivudine and zalcitabine showed promising inhibitory activity (Table S3). Finally, but very importantly, of the 1PPs, the chemotherapy drug uracil mustard showed binding to the catalytic dyad with all of its simulated conformations, followed by anti-cytomegalovirus cidofovir (Table S4). glaucoma citicoline, and anti-HIV drugs lamivudine and zalcitabine showed promising inhibitory activity (Table S3).
Finally, but very importantly, of the 1PPs, the chemotherapy drug uracil mustard showed binding to the catalytic dyad with all of its simulated conformations, followed by anti-cytomegalovirus cidofovir (Table S4). (a)

Molecular Dynamics Simulations
MD was performed on the hole complexes of the top 11 PPs candidates. Based on the conformer score energy from docking, the complex with the conformer with the lowest value and interacting with the 3CL pro active site was selected. The RMSDs were computed along the trajectories using the initial structure as a reference. Figure 3 shows that the binding of PPs significantly affected the equilibration states of 3CL pro , as the majority of the tested systems reached their equilibrium at around 100 ns. The PP-3CLpro complexes revealed relatively lower average values for the RMSDs, between 0.41 and 0.52 Å, throughout the simulation, clarifying their good behavior in forming stable complexes. Moreover, the fluctuations in the 3CL pro backbone residues were analyzed by means of the RMSF (Figure 4). The 3CL pro /PPs exhibited lower fluctuations, particularly at the active site. The fluctuations at the catalytic dyad and GLU 166 were minor, demonstrating the loss of flexibility at these regions upon binding to the PPs. Table 3 shows the superior stability of the PP-3CL pro complexes formed throughout the production runs; these results support the use of these PPs as 3CL pro inhibitors. Table 3. The binding interactions of the potential pyrimidonic pharmaceuticals/3-chymotrypsin-like protease 3CL pro complexes at different times throughout the production runs.

Uracil mustard
Molecules 2021, 26, x FOR PEER REVIEW 10 of 18 MD was performed on the hole complexes of the top 11 PPs candidates. Based on the conformer score energy from docking, the complex with the conformer with the lowest value and interacting with the 3CL pro active site was selected. The RMSDs were computed along the trajectories using the initial structure as a reference. Figure 3 shows that the binding of PPs significantly affected the equilibration states of 3CL pro , as the majority of the tested systems reached their equilibrium at around 100 ns. The PP-3CLpro complexes revealed relatively lower average values for the RMSDs, between 0.41 and 0.52 Å, throughout the simulation, clarifying their good behavior in forming stable complexes. Moreover, the fluctuations in the 3CL pro backbone residues were analyzed by means of the RMSF (Figure 4). The 3CL pro /PPs exhibited lower fluctuations, particularly at the active site. The fluctuations at the catalytic dyad and GLU 166 were minor, demonstrating the loss of flexibility at these regions upon binding to the PPs. Table 3 shows the superior stability of the PP-3CL pro complexes formed throughout the production runs; these results support the use of these PPs as 3CL pro inhibitors. MD was performed on the hole complexes of the top 11 PPs candidates. Based on the conformer score energy from docking, the complex with the conformer with the lowest value and interacting with the 3CL pro active site was selected. The RMSDs were computed along the trajectories using the initial structure as a reference. Figure 3 shows that the binding of PPs significantly affected the equilibration states of 3CL pro , as the majority of the tested systems reached their equilibrium at around 100 ns. The PP-3CLpro complexes revealed relatively lower average values for the RMSDs, between 0.41 and 0.52 Å, throughout the simulation, clarifying their good behavior in forming stable complexes. Moreover, the fluctuations in the 3CL pro backbone residues were analyzed by means of the RMSF (Figure 4). The 3CL pro /PPs exhibited lower fluctuations, particularly at the active site. The fluctuations at the catalytic dyad and GLU 166 were minor, demonstrating the loss of flexibility at these regions upon binding to the PPs. Table 3 shows the superior stability of the PP-3CL pro complexes formed throughout the production runs; these results support the use of these PPs as 3CL pro inhibitors. MD was performed on the hole complexes of the top 11 PPs candidates. Based on the conformer score energy from docking, the complex with the conformer with the lowest value and interacting with the 3CL pro active site was selected. The RMSDs were computed along the trajectories using the initial structure as a reference. Figure 3 shows that the binding of PPs significantly affected the equilibration states of 3CL pro , as the majority of the tested systems reached their equilibrium at around 100 ns. The PP-3CLpro complexes revealed relatively lower average values for the RMSDs, between 0.41 and 0.52 Å, throughout the simulation, clarifying their good behavior in forming stable complexes. Moreover, the fluctuations in the 3CL pro backbone residues were analyzed by means of the RMSF (Figure 4). The 3CL pro /PPs exhibited lower fluctuations, particularly at the active site. The fluctuations at the catalytic dyad and GLU 166 were minor, demonstrating the loss of flexibility at these regions upon binding to the PPs. Table 3 shows the superior stability of the PP-3CL pro complexes formed throughout the production runs; these results support the use of these PPs as 3CL pro inhibitors. MD was performed on the hole complexes of the top 11 PPs candidates. Based on the conformer score energy from docking, the complex with the conformer with the lowest value and interacting with the 3CL pro active site was selected. The RMSDs were computed along the trajectories using the initial structure as a reference. Figure 3 shows that the binding of PPs significantly affected the equilibration states of 3CL pro , as the majority of the tested systems reached their equilibrium at around 100 ns. The PP-3CLpro complexes revealed relatively lower average values for the RMSDs, between 0.41 and 0.52 Å, throughout the simulation, clarifying their good behavior in forming stable complexes. Moreover, the fluctuations in the 3CL pro backbone residues were analyzed by means of the RMSF (Figure 4). The 3CL pro /PPs exhibited lower fluctuations, particularly at the active site. The fluctuations at the catalytic dyad and GLU 166 were minor, demonstrating the loss of flexibility at these regions upon binding to the PPs. Table 3 shows the superior stability of the PP-3CL pro complexes formed throughout the production runs; these results support the use of these PPs as 3CL pro inhibitors. MD was performed on the hole complexes of the top 11 PPs candidates. Based on the conformer score energy from docking, the complex with the conformer with the lowest value and interacting with the 3CL pro active site was selected. The RMSDs were computed along the trajectories using the initial structure as a reference. Figure 3 shows that the binding of PPs significantly affected the equilibration states of 3CL pro , as the majority of the tested systems reached their equilibrium at around 100 ns. The PP-3CLpro complexes revealed relatively lower average values for the RMSDs, between 0.41 and 0.52 Å, throughout the simulation, clarifying their good behavior in forming stable complexes. Moreover, the fluctuations in the 3CL pro backbone residues were analyzed by means of the RMSF (Figure 4). The 3CL pro /PPs exhibited lower fluctuations, particularly at the active site. The fluctuations at the catalytic dyad and GLU 166 were minor, demonstrating the loss of flexibility at these regions upon binding to the PPs. Table 3 shows the superior stability of the PP-3CL pro complexes formed throughout the production runs; these results support the use of these PPs as 3CL pro inhibitors. MD was performed on the hole complexes of the top 11 PPs candidates. Based on the conformer score energy from docking, the complex with the conformer with the lowest value and interacting with the 3CL pro active site was selected. The RMSDs were computed along the trajectories using the initial structure as a reference. Figure 3 shows that the binding of PPs significantly affected the equilibration states of 3CL pro , as the majority of the tested systems reached their equilibrium at around 100 ns. The PP-3CLpro complexes revealed relatively lower average values for the RMSDs, between 0.41 and 0.52 Å, throughout the simulation, clarifying their good behavior in forming stable complexes. Moreover, the fluctuations in the 3CL pro backbone residues were analyzed by means of the RMSF (Figure 4). The 3CL pro /PPs exhibited lower fluctuations, particularly at the active site. The fluctuations at the catalytic dyad and GLU 166 were minor, demonstrating the loss of flexibility at these regions upon binding to the PPs. Table 3 shows the superior stability of the PP-3CL pro complexes formed throughout the production runs; these results support the use of these PPs as 3CL pro inhibitors. MD was performed on the hole complexes of the top 11 PPs candidates. Based on the conformer score energy from docking, the complex with the conformer with the lowest value and interacting with the 3CL pro active site was selected. The RMSDs were computed along the trajectories using the initial structure as a reference. Figure 3 shows that the binding of PPs significantly affected the equilibration states of 3CL pro , as the majority of the tested systems reached their equilibrium at around 100 ns. The PP-3CLpro complexes revealed relatively lower average values for the RMSDs, between 0.41 and 0.52 Å, throughout the simulation, clarifying their good behavior in forming stable complexes. Moreover, the fluctuations in the 3CL pro backbone residues were analyzed by means of the RMSF (Figure 4). The 3CL pro /PPs exhibited lower fluctuations, particularly at the active site. The fluctuations at the catalytic dyad and GLU 166 were minor, demonstrating the loss of flexibility at these regions upon binding to the PPs. Table 3 shows the superior stability of the PP-3CL pro complexes formed throughout the production runs; these results support the use of these PPs as 3CL pro inhibitors. MD was performed on the hole complexes of the top 11 PPs candidates. Based on the conformer score energy from docking, the complex with the conformer with the lowest value and interacting with the 3CL pro active site was selected. The RMSDs were computed along the trajectories using the initial structure as a reference. Figure 3 shows that the binding of PPs significantly affected the equilibration states of 3CL pro , as the majority of the tested systems reached their equilibrium at around 100 ns. The PP-3CLpro complexes revealed relatively lower average values for the RMSDs, between 0.41 and 0.52 Å, throughout the simulation, clarifying their good behavior in forming stable complexes. Moreover, the fluctuations in the 3CL pro backbone residues were analyzed by means of the RMSF (Figure 4). The 3CL pro /PPs exhibited lower fluctuations, particularly at the active site. The fluctuations at the catalytic dyad and GLU 166 were minor, demonstrating the loss of flexibility at these regions upon binding to the PPs. Table 3 shows the superior stability of the PP-3CL pro complexes formed throughout the production runs; these results support the use of these PPs as 3CL pro inhibitors. MD was performed on the hole complexes of the top 11 PPs candidates. Based on the conformer score energy from docking, the complex with the conformer with the lowest value and interacting with the 3CL pro active site was selected. The RMSDs were computed along the trajectories using the initial structure as a reference. Figure 3 shows that the binding of PPs significantly affected the equilibration states of 3CL pro , as the majority of the tested systems reached their equilibrium at around 100 ns. The PP-3CLpro complexes revealed relatively lower average values for the RMSDs, between 0.41 and 0.52 Å, throughout the simulation, clarifying their good behavior in forming stable complexes. Moreover, the fluctuations in the 3CL pro backbone residues were analyzed by means of the RMSF (Figure 4). The 3CL pro /PPs exhibited lower fluctuations, particularly at the active site. The fluctuations at the catalytic dyad and GLU 166 were minor, demonstrating the loss of flexibility at these regions upon binding to the PPs. Table 3 shows the superior stability of the PP-3CL pro complexes formed throughout the production runs; these results support the use of these PPs as 3CL pro inhibitors.

The Binding Free Energies
The data on the binding free energies of the PP-3CL pro complexes are tabulated in Table 4. In both MM-GBSA and MM-PBSA, van der Waals and electrostatic interactions

The Binding Free Energies
The data on the binding free energies of the PP-3CL pro complexes are tabulated in Table 4. In both MM-GBSA and MM-PBSA, van der Waals and electrostatic interactions acted as driving forces for the PP ligands to bind to 3CL pro , contrasting the solvation energies. The MM-GBSA and MM-PBSA results suggest that the PPs have an excellent ability to inhibit 3CL pro . Of these PPs, citicoline revealed the most promising inhibitory activity, followed by uridine triacetate (Table 4). Among the challenges of discovering 3CL pro inhibitors for COVID-19 treatment, these inhibitors must be highly bioavailable inside the cytosol [13]. The 2aPP and 2bPP structures contain a primidone heterocycle and ribose ring. These heterocycles increase their hydrophilicity and solubility in plasma. Thus, they can satisfy the requirement of bioavailability. Citicoline, the most promising inhibitor among the PPs investigated, has high hydrophilicity and good ADME properties [42,43]. The 2aPPs and 2bPPs also have an intermediate structure size among the PP groups. This sheds light on the importance of the size and general structural features of PPs acting as 3CL pro inhibitors.
Recent reports suggest a general hypothesis that 3CL pro inhibitors comprise electrophilic sites such as Michael acceptors [12]. That the pyrimidone ring is highly electrondeficient clarifies and confirms this hypothesis. The pyrimidone ring plays an essential role in PPs' inhibitory activity and has a high tendency to form hydrogen bonds, particularly with the GLU 166 residue. This may preclude the formation of the S1 pocket. Contacts between the pyrimidone ring and HIS 41 were observed in floxuridine, stavudine, and telbivudine. In all these cases, HIS 41 interacts with the oxygen of the pyrimidone group. Further, the electrophilic carbon, nitrogen, and oxygen in the pyrimidone ring were attracted to bind with the sulfur of the CYS 145 residue.
To conclude, the inhibitory effect on 3CL pro by PPs was investigated based on their ability to form hydrogen bonds and van der Waals interactions with the 3CL pro active side through molecular docking, MD simulations, and the calculation of binding free energy. The overall analysis revealed 11 candidates from the initial set of 42 investigated PPs are promising 3CL pro inhibitors. These include citicoline and uridine triacetate as the best choices, followed by telbivudine, trifluridine, lamivudine, cytarabine, stavudine, zalcitabine, tipiracil, floxuridine, and flavin mononucleotide. The interactions of PPs with the catalytic dyad and the active sides of 3CL pro of SARS-CoV-2 were comprehensively and thoroughly investigated. The pyrimidone ring was found to play an essential role in the PPs' inhibitory activity.