In Silico Prediction of Inhibitory Potential of a Punicalagin β-anomer Against SARS-COV-2 Main Protease (Mpro)

The pandemic caused by the new coronavirus has resulted in a global health emergency and has prompted an urgent need for new treatment strategies. No target-specic drugs are currently available for SARS-CoV-2, but new drug candidates targeting the viral replication cycle are being explored. A prime target of drug-discovery efforts is the SARS-CoV-2 main protease (Mpro). In this work, we identied a potential inhibitor for SARS-CoV-2 main protease using in silico methodologies. Molecular docking and molecular dynamics studies were carried out to ascertain the inhibitory action of α and β anomers of Punicalagin from fruit peel of Punica granatum against the Mpro protease. The molecular dynamics results revealed that the β-anomeric conguration of punicalagin allowed access to more hydrogen bonds and hydrophobic interaction leading to higher selectivity and specicity of β-anomer than α-anomer. Therefore, the β-anomer of Punicalagin could act as potential inhibitor against the main protease of SARS-CoV-2 and may act as a potential drug candidate.


Introduction
The outbreak caused by a new species of coronavirus led the World Health Organization (WHO) in early 2020 to declare the disease caused by this virus as a Public Health Emergency of International Importance and soon after it was characterized as a pandemic [1]. In December 2019, in the city of Wuhan, China, an outbreak occurred, in which several individuals went to local hospitals with severe pneumonia of unknown etiology. After investigations to nd out about the cause of this outbreak, it was found that most of the initial cases had contact with the same wholesale seafood market [2]. The etiologic agent discovered was a new virus of the family Coronaviridae and genus β-coronavirus and was identi ed as severe acute respiratory syndrome of coronavirus 2 (SARS-CoV-2), which causes coronavirus 2019 or COVID-19, as it became known. In addition, it was shown that, unlike other coronaviruses that infect humans, SARS-CoV-2 has a high human-to-human transmission rate, which was a determining factor for its dissemination [3,4].
The virus can infect individuals of all ages, however it mainly affects men over 60 years [5]. According to He,Deng and Li[6], the lethality rate is higher in patients with coexisting medical conditions, such as diabetes, hypertension and cardiovascular diseases. Symptoms vary from fever, cough, sore throat, changes in smell and taste to acute respiratory distress syndrome (ARDS) and dyspnoea in the most severe cases[6-8].
The SARS-CoV-2 virus consists of four structures, the spike glycoprotein or S glycoprotein, a membrane protein, an envelope protein and a nucleocapsid protein, which is present within the viral envelope [9].
SARS-CoV-2 is a ribonucleic acid (RNA) virus whose genetic material is a single positive RNA molecule. The virus enters the human host, through the binding of the spike glycoprotein, in the S1 subunit, with its receptor, the angiotensin-converting enzyme 2 (ACE2) that is present on the cell surface [10]. After binding, Furin enzyme will pre-cleave the S glycoprotein in the S1 and S2 subunits and then the transmembrane serine protease 2 (TMPRSS2) enzyme will cleave the S glycoprotein in the S2 subunit. Subsequently, the S2 subunit will assist in the fusion of the viral membrane with the host cell membrane and thus, with the entry of the virus into the host cell [4,11].
The SARS-CoV-2 main protease is an enzyme called Mpro or 3C-like protease (3CLpro), that plays an important role in the replication cycle, thus becoming a potential target for the control of viral infection [12]. Its action occurs after the virus enters the host cells, where the viral RNA genome will translate inactive proteins, called pp1a and pp1ab, which will be cleaved through Mpro and will result in non-structural proteins and join the RNA to form the virus genome [13]. Thus, it becomes necessary means for inhibiting replication of the SARS-CoV-2 virus. Researchers are seeking out products with e cient antiviral activity to prevent the SARS-CoV-2 virus replication. The peel of pomegranate (Punica granatum), has substantial studies that encompass antibacterial [14], anti-in ammatory [15] and antioxidant [16] activities in addition to antiviral activity agaist In uenza virus [17], Human Immunode ciency virus (HIV)[18] and Herpes simplex virus [19]. Punicalagin is a phenolic compound present abundantly in pomegranate peel with high molecular weight [20]. Punicalagin is found naturally in two anomeric forms α and β differing from the chiral carbon position (Fig. 1). Recent studies including α/ β anomers [21][22][23] revels that there are differences in the non-covalent interactions pro le between these anomers and proteins.
Hence, the present in silico study was aimed to explore the differences in inhibitory potential of α/β anomers of Punicalagin against the main protease viral of SARS-CoV-2 (Mpro).

Preparation of Mpro and α/β anomers
The protonation states of the molecules has a signi cant in uence on its mode of interaction with the receptor. Based on the pKa values, the protonation states of the α/β anomers studied were calculated using the ChemAxon's MARVIN software [24].
The geometries of the α/β anomers were optimized by using standard techniques [25]. The X ray crystal structure of COVID-19 main protease (Mpro) in complex with an inhibitor N3 [31] at 2.16 Å resolution (PDB ID: 6LU7) used is this study was obtained from RCSB Protein Data Bank [32]. The Mpro structure was prepared for molecular docking by removing the water molecules and other hetero molecules from the original crystal structure.
The molecular docking study Docking Input les were created using AutoDock tools of MGL tools [33,34]. The 3D grid box with a box size of 30 Å x 30 Å x 30 Å with the center coordinates x = -13.579, y = 18.940 and z = 69.318, which are the coordinates in the Nε2 of His41 residue that is included in the Mpro catalytic dyad [35]. Molecular docking was performed using Autodock Vina [36]. Nine best poses were generated for each anomer and scored using Autodock Vina scoring functions. The lowest energy pose of each anomer was used as input for molecular dynamics (DM) simulations.

Molecular dynamics simulations
All molecular dynamics (MD) simulations were performed using the GROMACS 2018.4 [37] package implemented with the AMBER ff99SB[38]. The transferable intramolecular potential with 3 points (TIP3P) [39] water molecules were used to solvate the simulated systems. The systems neutralization was achieved through the addition of counter ions. The Leap-Frog algorithm [40] was applied to integrate the motion equation with time step of 2.0 fs. The long-range interactions were modeled using particle-mesh Ewald sum (PME) [41] with a cut-off of 1.2 nm. The van der Waals interactions were also calculated using the same threshold. Bonds involving hydrogen atoms were restrained using LINCS algorithm [42]. The Nosé-Hoover thermostat [43] was used to x the system temperature (310 K) in all production simulations, while the system pressure was controlled using a Parrinello-Rahman barostat [44] in the NPT simulations. The geometry of the systems were minimized by the steepest descent algorithm [45] for 10000 steps with tolerance of 10 kJ mol − 1 nm − 1 followed by conjugate gradient algorithm[46] for 10000 steps with tolerance of 10 kJ mol − 1 nm − 1 . Two shorts 20-ns equilibrium dynamics with NVT and NPT ensembles were performed. Finally, 200-ns production MD simulation using NVT ensemble was performed in three replicates for each system in order to determine the interaction energy between α/β anomers and Mpro.
To analyze the interaction energy between α/β anomers and Mpro we used the Interaction Potential Energy (IPE) [47], which can be de ned as the total interaction energy between two groups, α/β anomers and Mpro. The IPE calculation includes the sum of van der Waals and electrostatic contributions. To identify non-covalent interactions between anomers and Mpro receptor, we used an open-source proteinligand interaction pro ler (PLIP) [48]. The graphical representations of the protein-receptor interaction were built with Visual Molecular Dynamics (VMD) [49] software from PLIP analysis.

Results And Discussion
The pKa can predict the protonation states of organic molecules across pH range. Ionizable sites are found often in organic molecules and in uence their pharmaceutical properties including target a nity [50]. Due to the importance of knowing the protonation state, the microspecies distribution (%) as a function of pH (1.0-14.0) for α/β anomers were performed (Supplementary Fig. S1). The pKa values of hydroxyl groups are shown in the Supplementary Figs. S2a (α anomer) and S2b (β anomer). We can see that the pKa values do not change with con guration change and then α/β anomers have the same ionizable sites at pH = 7.4, that are indicated by blue arrows in Supplementary Fig. S2a. Hence, negatively charged punicalagin α/β anomers were considered for this study. In order to improve the structures of α/β anomers their optimized structures were performed by Density Functional Theory (DFT) using B3LYP/6-31 + G(d,p) level of theory (Fig. 2).
From α and β optimized structures and using Mpro as receptor, the best poses (pose 1) generated from docking analysis (Supplementary Fig. S3), MD simulations were performed for identifying the α/β anomers intermolecular interaction pro le. The interaction potential energy (IPE) calculations were performed to evaluate the a nity as well as the contribution of each residue of binding pocket of Mpro receptor with the ligands. Based on atoms of α/β anomers as reference, the RMSD analysis (Fig. 3) showed a structural equilibrium around 125 ns for α (Fig. 3A, black, red and green lines) and around 50 ns for β (Fig. 3B, blue, orange and purple lines). Thus, the IPE calculations were performed in the last 75 ns for α anomer and 150 ns for the β anomer. The average IPE (standard deviation) between Mpro receptor and α/β anomers were − 250.96 kJ mol − 1 (20.40) and − 320.00 kJ mol − 1 (27.05), respectively. Then, according to the IPE values, β anomer has a highest a nity for Mpro compared to the α anomer. This result suggests that the β anomer can become a potential inhibitor for SARS-CoV-2 main protease.
A closed comparative analysis of the binding of α (Fig. 4a) and β (Fig. 4b) anomers with pocket Mpro residues in the last simulation frame (200 ns) of MD, revealed that the β anomer binds to the protein through nine favorable hydrogen bonds. On the other hand, a close inspection of the binding of α anomer ( Fig. 4a) with pocket Mpro residues showed that α anomer ligand binds with the protein with seven hydrogen bonds. Hydrogen bonds contribute from 11 to 60 kJ mol − 1 in the interaction potential energy (IPE) [51] and are force intermolecular most in uential in molecular recognition [52]. Furthermore, the βanomeric carbon con guration (Fig. 5) allowed the holding of two hydrogen bonds (blue dashed line) and one hydrophobic interaction (gray dashed line) with THR45 Mpro residue. We can see in Fig. 5 that two hydrogen bonds are formed between hydroxyl groups of β-anomeric carbon and residue THR45 at distances of 3.42 Å and 4.74 Å. In addition, the hydrophobic interaction is formed between C-γ of residue THR45 and C-47 of β anomer. The presence of hydrophobic groups in the ligand can expel molecular water from protein pocket and maximize hydrogen bonds account leading to an increased selectivity and e ciency of the receptor for a certain ligand [53], as in the case of the β-anomer and Mpro receptor, where the β-anomeric con guration of punicalagin allowed access to two more hydrogen bonds and a hydrophobic interaction. In other words, the combination of hydrogen bonds and hydrophobic interaction found between β-anomer and Mpro binding pocket led to higher selectivity and speci city of the βanomer ligand compared to the α-anomer.

Conclusion
In this work, we reported a systematic study of a nity and speci city of α-and β-anomers of Punicalagin with Mpro protease of SARS-CoV-2 by using molecular docking and molecular dynamics simulations.
Interaction Potential Energy (IPE) calculations through molecular dynamics simulations reveals that the β-anomer of Punicalagin has a higher a nity for the pocket Mpro receptor compared to the α-anomer. Furthermore, β-anomer con guration shows more hydrogen bonds and a hydrophobic interaction than αanomer with pocket Mpro residues. The amount and combination of hydrogen bonds and hydrophobic interaction can lead to selectivity and speci city of β-anomer of Punicalagin being able to act as a potential drug candidate. Therefore, we way conclude that the β-anomer ligand could act as potential inhibitor against the main protease of SARS-CoV-2. Code availability Not applicable.

Data Availability
All data are available on request to the corresponding author.

Con ict of interest
The authors declares no competing interests. The chemical structures of α and β anomers of Punicalagin.

Figure 3
Root Mean Square Deviations (RMSD) for punicalagin α (A, black, red and green lines) and β (B, blue, orange and violet lines) anomers. All MD simulations were performed in three replicates each anomer studied.