Mechanistic Insights into the Inhibition of SARS-CoV-2 Main Protease by Clovamide and Its Derivatives: In Silico Studies

: The novel coronavirus SARS-CoV-2 Main Protease (M pro ) is an internally encoded enzyme that hydrolyzes the translated polyproteins at designated sites. The protease directly mediates viral replication processes; hence, a promising target for drug design. Plant-based natural products, especially polyphenols and phenolic compounds, provide the scaffold for many effective antiviral medications, and have recently been shown to be able to inhibit M pro of SARS-CoV-2. Speciﬁcally, polyphenolic compounds found in cacao and chocolate products have been shown by recent experimental studies to have strong inhibitory effects against M pro activities. This work aims to uncover the inhibition processes of M pro by a natural phenolic compound found in cacao and chocolate products, clovamide . Clovamide (caffeoyl-DOPA) is a naturally occurring caffeoyl conjugate that is found in the phenolic fraction of Theobroma Cacao L. and a potent radical-scavenging antioxidant as suggested by previous studies of our group. Here, we propose inhibitory mechanisms by which clovamide may act as a M pro inhibitor as it becomes oxidized by scavenging reactive oxygen species (ROS) in the body, or becomes oxidized as a result of enzymatic browning. We use molecular docking, annealing-based molecular dynamics, and Density Functional Theory (DFT) calculations to study the interactions between clovamide with its derivatives and M pro catalytic and allosteric sites. Our molecular modelling studies provide mechanistic insights of clovamide inhibition of M pro , and indicate that clovamide may be a promising candidate as a drug lead molecule for COVID-19 treatments. a native we that clovamide could the a transfer We suggest by two of the is also capable of inhibiting by inducing local enzymatic and altering hence, inhibiting activity by inhibition We hope this work future experimental studies to test our computational this a future structure-based rational design of protease inhibitors and future developments of


Introduction
The main protease (M pro ) of SARS-CoV-2, also known as the Chymotrypsin-like cysteine protease (3CL pro ), is the critical protease of the coronavirus and bears the essential role of mediating viral replication and transcription processes, making it one of the most attractive drug targets for SARS-CoV-2 [1,2]. Coronaviruses are subject to extensive mutagenesis; however, key proteins such as the main protease are highly conserved, since mutations in these key proteins are often lethal to the virus. For this reason, drugs that target the conserved M pro are often able to prevent viral replication and proliferation processes, and display a broad-spectrum antiviral activity [3].
A large number of crystal structures of M pro from SARS-CoV-2 has been solved, either in their apo forms, or in complex with different enzyme inhibitors, providing important information pertaining protein folding, assembly, and the catalytic mechanisms [4,5]. In particular, the SARS-CoV-2 M pro -N3 inhibitor complex was initially determined, suggesting a Michael addition mechanism that prompts the formation of the covalent bond between the deprotonated C145 residue and the α,βunsaturated carbonyl moieties [6]. The N3 inhibitor was further modified with α-ketoamide functional groups to improve both the efficacy of the inhibitor and its movement across the cell membrane [7].
Many effective pharmaceuticals are derived from plant-based natural compounds. These include antiviral drugs for the treatment of various diseases by inhibiting the virus replications, such as HIV [8] or Hepatitis C [9]. Recent studies revealed a wide range of natural compounds that can serve as antiviral drugs for COVID-19 treatments; in particular, many reports have detailed natural phytochemicals that act to inhibit SARS-CoV-2 M pro activities [10][11][12]. Polyphenols and phenolic compounds are a major class of bioactive natural compounds and are known for their antiviral and pleiotropic effects. Many in silico and in vitro studies have suggested that phenolic compounds found in plant metabolites are capable of inhibiting SARS-CoV-2 M pro , such as resveratrol, myricetin, quercetin, amentoflavone, and puerarin [13][14][15][16][17]. Recent work by Zhu and Xie (2020) showed that cocoa and dark chocolate extracts are capable of inhibiting the activity of SARS-CoV-2 M pro [18].
Clovamide (trans-clovamide, (2S)-3-(3,4-dihydroxyphenyl)-2-[[(E)-3-(3,4-dihydroxyp henyl)prop-2-enoyl]amino]propanoic acid) is a naturally occurring caffeoyl conjugate ( Figure 1) initially identified in red clover, from which its name was derived, and later found in Theobroma Cacao L. [19]. Though only found in a minor quantity in cacao phenolic fractions, clovamide is a potent antioxidant with neuroprotective effects against lipid oxidation and reactive oxygen species (ROS) generation [20]. As such, it is a potential candidate for the therapeutic treatment of Parkinson's disease [21]. Previous work by our research group identified and quantified clovamide in cocoa beans at different production stages, and has shown the impact of high-temperature processing on clovamide content; we also employed an electrochemical method and confirmed its outstanding antioxidant property [22]. This present work explores the potential inhibitory and antiviral effects of clovamide on SARS-CoV-2 M pro . Previously, clovamide was reported to be a potent inhibitor of proteolysis in forage crops, as a result of its oxidation and crosslinking with proteins [23]. We hypothesized that clovamide could act as a strong inhibitor to the main protease, for the following reasons: (1) Two catechol (o-hydroxy phenol) moieties present in the clovamide molecule, which are known to react with cysteine residues readily in oxidized states; when clovamide scavenges radical species, the catechol groups become oxidized into orthoquinones [24,25]. Alternatively, the oxidation and inhibition can be mediated by natural enzymatic browning or melanization processes, generating quinone species by polyphenol oxidase (PPO)-mediated oxidation and subsequent crosslinking with proteins [26,27]. The non-enzymatic oxidation of o-diphenol to quinone species of clovamide was also detected [26]. Quinones, particularly ortho-quinones, are shown to be reactive to thiolcontaining compounds which gives them potent antiviral activities [28,29]. (2) There is a α,βunsaturated carbonyl moiety. Our previous study determined its essential role in providing conjugation that stabilizes the radical scavenging process by one of the two catechol moieties [22]. This α,βunsaturated carbonyl may serve as the acceptor in the Michael addition reaction as described in the N3 inhibitory mechanism [6]. (3) Other essential functional groups in the molecule could facilitate its binding in the SARS-CoV-2 M pro substrate-binding pocket, such as the carboxylic acid group, via non-covalent interactions such as hydrogen bonding [30].
The antiviral effect of clovamide was first reported by El-Sharawy et al. (2017) against the H5N1 influenza A virus [31]. Recently, clovamide was suggested to be a broad-spectrum enzyme inhibitor adopting a generic mechanism involving the o-diphenol oxidation to its quinone form and indiscriminate covalent cross-linking to reducing groups, such as thiols or amines, in proteins [26].
Developing a lead molecule and an effective drug (small molecules with desired medicinal properties) is challenging even with known biological targets [32]. At the early stage of drug development, rational drug design is an essential tool for the preclinical evaluations of drug properties of compounds [32]. Modern computing technology enables the use of chemo/bioinformatics and molecular modelling tools to study the binding activities of drugs and medicines with desired targets. Molecular docking can identify lead molecules against target enzymes [33,34], and molecular dynamics (MD) aid the study of protein-ligand interactions by exploring the binding affinities under different energetic conditions. Density Functional Theory (DFT) and quantum mechanical studies closely examine the enzyme catalytic processes at the active site and the mechanisms by which the inhibitor prevents the catalytic activity of the enzyme [35]. Together, they provide us with mechanistic insights regarding how small molecules could act as enzyme inhibitors.
This work explores the potential inhibitory and antiviral effects of clovamide and its structural analogues, on SARS-CoV-2 M pro . We used in silico techniques introduced above to investigate the binding affinities and potential mechanisms based on which these small molecules could inhibit the main protease catalytic activities and, therefore, stop the viral replication process. Our work suggests that clovamide is a potential main protease inhibitor and we propose conjugate addition mechanisms by which clovamide analogues, upon radical scavenging or natural enzymatic browning, could act as irreversible inhibitors to inhibit M pro catalytic activities. We also demonstrate that clovamide could act as an allosteric reversible inhibitor, inhibiting the enzyme by inducing protein allostery. We hereby show that clovamide may be a promising candidate as a drug lead molecule for COVID-19 treatments by targeting the critical M pro and inhibiting viral replication processes.

Density Functional Theory
Density Functional Theory (DFT) was applied for energetic and geometric calculations of different reaction mechanisms to gain thermodynamic insights. We performed these quantum mechanical calculations using Unix-based Q-Chem (version 5.3.0) program [36]. Molecular models were built using IQmol (version 2.15.0) and was used to generate input files for Q-Chem. The hybrid functional B3LYP and basis set 6-31 G were used. Self-consistent field (SCF) convergence tolerance was set to 1 × 10 −7 . Dispersion Correction DFT-D3(BJ) was included in the calculations. An implicit solvent model PCM was employed with water (dielectric constant = 78.39). This combination of DFT-D method and PCM model was shown to be very efficient for accurate treatments of non-covalent interactions in large, solvated, non-covalently bound environments, especially protein and other biologically relevant systems [37]. In our calculations, we adopted the QM-cluster approach [35,38], and selected the ligand (the inhibitor clovamide) with its surrounding amino acids in the enzyme active site that were within 5.0 Å away from the center of the inhibitor molecule. This included the following 22 amino acids, accounting for a total of 380 atoms: T25, T26, L27, H41, V42, M49, L141, N142, G143, S144, C145, G146, H163, H164, M165, E166, L167, P168, Q189, T190, A191, and Q192. To preserve the protein environment while providing chemical accuracy, we selectively constrained the 66 peptide backbone atoms (α-carbons, backbone carbonyl carbons, and backbone amine nitrogen atoms of all residues). All other atoms (side chain atoms) of amino acid residues and the entire ligand molecule were allowed flexible. The ends of these short peptide segments were capped with amine (-NH 2 ) groups at N-terminus and hydroxy (-OH) groups at C-terminus to maintain neutral charge and polarity found in peptide-bonded residues in the full enzyme [39].

Molecular Docking and Dynamics
Molecular docking and dynamics calculations were performed using software from BIOVIA (San Diego, CA, USA) with CDOCKER package in Discovery Studio 2020 version [40]. For molecular dynamics, we used default conditions set in the Standard Dynamics Cascade protocol.
The SARS-CoV-2 M pro was obtained from the PDB database (PDB ID: 6LU7) with the inhibitor N3 co-crystallized. Using BIOVIA Discovery Studio 2020, we followed the standard procedure of preparing this protein that provides assignment of the CHARMM force field and the addition of hydrogen atoms. The inhibitor N3 was removed, and the binding site was chosen after selecting the catalytic dyad amino acids, for a sphere of radius 18 Å. Molecular docking was performed for 15 poses. Then, the docking pose that was of our interest (ad hoc for each specific docking study) was selected for Molecular Dynamics (MD) study. We applied an integrated Standard Dynamics Cascade protocol incorporated in the Discovery Studio program, which includes five steps: (i) solvation, (ii) energy optimization (minimization), (iii) built-in annealing process (heating and cooling), (iv) equilibrium run, and (v) production run. The MD results provided us with closer examinations of how the ligand binding and inhibitory activities evolve in non-rigid protein environments. Final outputs from DFT, molecular docking, and molecular dynamics calculations were visualized and analyzed using Discovery Studio Visualizer and PyMOL (version 2.4.0, Schrödinger, LLC, New York, NY, USA).

SARS-CoV-2 M pro and Clovamide Molecule
Based on our hypotheses of the potential inhibitory mechanisms (Scheme 1), we started by exploring the binding activities of the clovamide molecule to the M pro active site. Results of docking were analyzed focusing on the interactions between S(C145) and the centroid one of the catechol rings of the clovamide molecule, which are found to be closely related in pose 1 and 12. Pose one had a CDOCKER interaction energy of −45.8 Kcal/mol and was selected for further study: the Discovery Studio molecular mechanics program provides the CDOCKER interaction energy, which is the non-bonded interaction energy (composed of the van der Waals term and the electrostatic term) between the protein and the ligand related to the force field CHARMM [40]. Figure S1 (Supplemental Information) shows essential amino acids involved in this interaction. The S(C145) was close to the centroid of the catechol ring, with a distance of 2.987 Å, providing a strong S-π interaction that stabilized the protease-inhibitor complex. To further explore the potential reaction mechanism, this docking pose was treated with a standard dynamic cascade and in situ ligand minimalization protocols. The dynamic cascade provided no promising results, in which the C145 residue evolved to be far apart from both catechol moieties. Figure 2 illustrates the result of the in situ minimalized interactions, with an increased S(C145) to centroid distance of 5.330 Å. However, instead of the anticipated Michael addition mechanism, we found a potential inhibition mechanism by which the deprotonated cysteine residue attacked the carbonyl carbon of the carboxylic acid functional group C(COOH), forming a thioester group. The minimized structure showed O(H163) hydrogen bonds to H-S(C145) of 2.795 Å, suggesting hydrogen capture by the oxygen and the thiolate formation, for further attack on the clovamide C(COOH) by the latter. The S(C145) was in proximity to the C(COOH) of clovamide, 3.110 Å. In fact, this motif of cysteine attacking the C(COOH) of clovamide molecule forming thioester was quite common in our computational results, which are described in a later section. Finally, O(F140) formed hydrogen bonds to both of the hydroxy groups of the catechol ring to stabilize the complex, with distances between O(F140) to the two hydroxy groups 2.517 Å and 2.788 Å, respectively. Scheme 1. Proposed inhibitory mechanisms based on Michael addition conjugation reactions. The cysteine residue can be deprotonated either by the later negatively charged oxygen (mechanism (A)), or by the H41 side-chain imidazole group (mechanism (B)) [6]. The critical C145 sulfur atom is colored in blue and the hydrogen atom in red to indicate their fate throughout the reactions. Regardless of the inhibitory pathways, the same covalently bound adduct forms. A Michael addition-based inhibitory mechanism we hypothesized earlier is illustrated in Scheme 1. However, our docking and dynamics studies did not show evidence suggesting this reaction to occur. Docking poses did not show proximity of the C145 residue to the double-bond β carbon.

SARS-CoV-2 M pro and Ortho-Quinone Derivatives of Clovamide
Clovamide is a natural phenolic antioxidant, whose potent antioxidant property is based on sequential radical reactions, including radical trapping from radical species generated from biomolecules. It has been shown that, under radical oxidation conditions, phenolic radicals can react with thiol radicals to produce sulfur-bearing phenolic compounds [41]. This modification is likely to happen under conditions that cause oxidative stress that produces ROS. Furthermore, Fujimoto and Masuda (2012) reported the forma-tions of thiol adducts of phenolic acid esters, flavonoids, and polyphenols when reacted with a cysteinyl thiol derivative, N-benzoylcysteine methyl ester, and DPPH to carry out radical oxidation reactions [42]. In this study, we hypothesized that under oxidative stress, phenolic radicals may form, and that catechol groups may undergo oxidation to form orthoquinone (Scheme 2) [42]. Alternatively, natural enzymatic browning can also produce quinone derivatives of clovamide [27].
The oxidized quinone species can then readily react with sulfur radicals or thiolate species of that thiol-containing cysteine residue [26]. Therefore, we investigated the three ortho-quinone derivatives ( Figure 3) of clovamide to explore their potential reactivities with the C145 residue in forming covalent bonds. A previous study by our group showed that between the two catechol rings in the clovamide molecule, the one that connected to the C=C double bond had a greater capability in scavenging superoxide anion radicals, due to the extra resonance provided by the conjugation [22]. We, therefore, proceeded to explore this clovamide analogue first ( Figure 3A).  Figure 4 shows the resulting conformation of the dynamic cascade of the clovamide quinone derivative A ( Figure 3A). In this conformation, the cysteine residue was close to the quinone ring, with a sulfur-centroid distance of 3.450 Å ( Figure 4A). We propose that this relative proximity could induce inhibitory pathways via different mechanisms. First, the deprotonation of the thiol hydrogen of the C145 can be achieved by backbone carbonyl oxygen of the H163 residue as shown previously [44]. In this conformation, the distance between H(C145) and O(H163) was 4.774 Å, rendering it unable to capture the hydrogen. However, by adjusting the C145 torsion angle C-C-S-H, the thiol hydrogen H(C145) moved closer to the ortho-quinone O(carbonyl). Thus, O(carbonyl) could capture H(C145) and form a hydrogen bond with the side chain H(S144, hydroxy), which was most likely to form strong hydrogen bonds with carbonyl oxygen, as seen by the shortest bond length (1.998 Å). Therefore, upon hydrogen capture, the protonated carbonyl group becomes a stronger electrophile, ready for nucleophilic attacks. This attack would then be performed by the negatively charged deprotonated thiolate group, attacking either the carbonyl carbon or the α carbon ( Figure 4). When attacking the carbonyl carbon, the formation of the S-C covalent bond leads to a tetrahedral product where the bond forms at the same carbon to which the hydroxy group is attached, by a mechanism described in a previous study on the inhibition of SARS-CoV-2 M pro by a methide quinone Celastrol [45]. Scheme 3(A1) illustrates the mechanisms. Given the challenges of hydrogen capture by O(H163), we proposed that the mechanism represented by the red arrow was favored.
When the thiolate ion attacked the α carbon, a Michael-type conjugate addition occurred, where the mono-thiol adduct formed as the product. The S-C covalent bond formed at the α carbon. Scheme 3 shows our proposed mechanism towards this reaction [42]. By the mechanism described in Scheme 3(A2), the major product where the addition occurred at the α carbon was preferred. This initial addition was, also, stereochemically governed; addition was less likely to happen at the carbon adjacent (and formed the minor product) to where the R group was located due to the resulting steric clash (Scheme 3B). The conjugate addition reaction that happens at the β carbon, mechanism of which is suggested in Scheme 3(A3)), is much less favored [29].
The entire inhibitor-protease complex was stabilized by numerous non-covalent interactions dominated by hydrogen bonds. Shown in Figure 4B,C is the 2D display of noncovalent interactions of the docking pose after dynamics cascade. Five pairs of hydrogen bonds stabilized the inhibitor-protease complex and the close proximity of C145 to the centroid of the quinone ring provided a π-π interaction. Three potential pairs of hydrogen bonds at the quinone carbonyl group made it a strong electrophile, to be attacked by the thiolate ion.
Next, we explored the clovamide quinone derivative B ( Figure 3B), shown in Figure 5. In this derivative, the ortho-quinone group replaced the catechol ring at the side without an extended conjugation. Results of the docking (15 poses) were analyzed focusing on interactions that could induce reactions described in Scheme 3. This conformation ( Figure 5A) of pose one showed a promising proximity between S(C145) and quinone centroid of 3.602 Å. The CDOCKER interaction energy for this conformation was −46.9 Kcal/mol. Figure 5A shows the docking results of pose 1. The docking pose showed that the H(C145) could be captured by backbone O(H163), as seen previously with the distance between H(C145) and O(H163) of 2.794 Å. In addition, adjusting the C145 torsion angle C-C-S-H from −60 • to 140 • suggested the potential hydrogen capture by the O(carbonyl) adjacent to the α1 carbon, having O(C=O, quinone) and H(C145) distance of 2.730 Å; or it could be captured by the O(carbonyl) adjacent to the α2 carbon, with a O(C=O) to H(C145) distance of 2.946 Å. In this situation, both thiolate formation mechanisms described in Scheme 3(A1) were favored. A close examination of the distances between S(C145), the carbonyl carbons, and the two α carbons revealed different possible inhibitory mechanisms. The fact that there were many potential inhibitory pathways suggested the high likelihood of clovamide inhibition. We proposed same mechanisms explained in Scheme 3(A1,A2) were potential pathways of inhibition by derivative B. Figure 5B is the 2D display of noncovalent interactions that stabilized the complex, with a substantial number of hydrogen bonds and van der Waal interactions acting as major stabilization forces. This docking pose was subsequently treated with a standard dynamics cascade protocol; however, the result of the dynamics cascade did not show expected interactions. The S(C145) to centroid distance became 4.821 Å and the O(H163) to H(C145) distance evolved to 6.078 Å. However, the non-covalent interactions of the structure after dynamics cascade suggested a stable conformation of the complex ( Figure 5C). conjugate addition mechanisms. (A1) Arrow-pushing mechanisms of the C145 thiolate attacks the carbonyl carbon, forming a tetrahedral product. The blue arrows represent the mechanism in which C145 is deprotonated by H163 backbone carbonyl oxygen [44,45], followed by the attack of S(C145) on the carbonyl of the ortho-quinone, whose O(C=O) is protonated by hydrogen capture of the surrounding Ser 144 residue. The red arrow indicates the mechanism in which the quinone carbonyl oxygen draws the hydrogen of the C145 residue, deprotonating it to become a thiolate which, subsequently, attacks the carbonyl carbon; (A2) the arrow-pushing mechanism of thiolate ion (deprotonation steps omitted for simplicity) adds to the α carbon. The scheme on the bottom shows that addition via the other α carbon is also feasible. (A3) The arrow-pushing mechanism of thiolate ion (deprotonation steps omitted for simplicity) adds to the β carbon. (B) The initial addition of thiolate to ortho-quinones to form mono-thiol adduct possesses regiospecificity, which is both mechanistically and stereochemically governed. Additions at the α carbons are found to be more likely than β carbon additions. Statistical data obtained from [29]. Finally, we examined clovamide quinone derivative C ( Figure 3C) in which both catechol rings were oxidized into ortho-quinone groups ( Figure 6). In a similar fashion, results of docking were analyzed seeking for insights of the inhibitory reaction mechanisms. The docking results of this "double-quinone" molecule showed two promising poses: pose 1 and pose 15. In pose 1, the S(C145) was in proximity to the carboxylic acid functional group in the substrate molecule, a motif observed previously in multiple settings. The CDOCKER interaction energy was −41.3 Kcal/mol. Pose 15 showed a promising ligand-protein complex structure in which the S(C145) pointed towards the center of the quinone ring at the conjugation side with an atom-to-centroid distance of 2.997 Å. Similar to derivative B, studies of this quinone analogue revealed that the deprotonation of the C145 residue could happen via hydrogen capture by either O(H163) or the carbonyl oxygen of the quinone group. When the C145 residue C-C-S-H torsion angle was −60 • , hydrogen was captured by backbone O(H163) with a distance between H(C145) and O(H163) of 2.792 Å; when the C-C-S-H torsion angle was adjusted to 105 • , hydrogen capture was likely to happen by O(C=O) of the carbonyl adjacent to carbon α2 with a distance of 2.756 Å, or by the oxygen of carbonyl group adjacent to carbon α1, with a distance between O(C=O) and H(C145) of 2.945 Å. The distances suggested that either proton capture mechanism could happen (described in Scheme 3(A1)). In addition, this docking pose also revealed multiple potential inhibitory mechanisms, including thiolate ion attacking either of the carbonyl carbons or either of the α carbons. The CDOCKER interaction energy for pose 15 was −35.1 Kcal/mol. Figure 6A shows the docking conformation of pose 15, which details all the distances that suggest reaction mechanisms. Figure 6B is the 2D display of non-covalent interactions that keep the inhibitor-protease complex stabilized.
To further investigate whether the conformations we obtained from molecular docking and dynamics studies were favored in real enzymatic environments, we performed DFT geometry and free energy calculations with the implicit solvent model and dispersion corrections (see Section 2.1). The results of geometry optimization confirmed that all three mechanisms proposed in Scheme 3 were thermodynamically favored, supporting our mechanistic proposal for clovamide inhibitions. The DFT free energy calculation suggested that the second α carbon mechanism (Scheme 3(A2)) was the most energetically stable product while the carbonyl carbon mechanism (Scheme 3(A1)) was the least stable one. All converged enzyme-inhibitor complexes showed relaxed conformations to accommodate evolved protein environments and a newly formed covalent S-C linkage ( Figure S2). Table S1 (Supplemental Information) lists the final energy of the starting structure and products we used for this evaluation.

Carboxylic Acid Is the Key Substrate Functional Group for C145 Thiolate Attacks
Our docking and dynamics results revealed a recurring structural motif in which the C145 residue was in proximity to the carboxylic acid functional group in the clovamide (or its derivative) molecules ( Figure S3). We proposed that carboxylic acid was a key functional group to which C145 thiolate ion could attack and form a thioester covalent linkage. This is a common reaction that can be found in many biological pathways; for instance, acyl groups are linked to coenzyme A (CoA) by thioester bonds, forming acyl CoA: an acyl group that is often linked to CoA is the acetyl unit, resulting in the derivative of acetyl CoA [46]. Figure S3 details the conformations and the distances between S(C145) to the C(COOH) of the clovamide (and derivative) compounds in docking studies.
To explore the potential inhibitory mechanism via the reaction between C145 and the carboxylic acid moiety of the clovamide molecule, we re-examined the canonical molecular docking outcome between clovamide and M pro . Pose 12 (orange compound shown in Figure S3) was selected for this study. Figure 7A shows the result of this docking study. As shown in previous cases, the deprotonation of the thiol hydrogen of C145 was expected to be completed by the backbone carbonyl oxygen of H163; the H (C145) to O(H163) distance was 2.796 Å. Following deprotonation, the resulting thiolate ion was, then, ready to attack the carbonyl carbon of the carboxylic acid group in a distance of 3.150 Å, and achieved inhibition by a thioesterification mechanism (proposed in Scheme 4). Initially, the carboxylate carbonyl oxygen became protonated by either the solvent environment, or the hydrogen of C145. When the C145 was at the conformation in Figure 7A where the C-C-S-H torsion angle was −60 • , the H(C145) to O(COOH) distance was 4.154 Å. When this torsion angle became 130 • , the distance was shortened to 1.760 Å, prone to be captured by the carboxylic acid carbonyl oxygen. We also observed a close-distance hydrogen bond between one imidazole nitrogen of H41 and the carboxylic acid carbonyl oxygen, 2.066 Å. Followed by the protonation of the acid were the nucleophilic attack by the thiolate ion, proton transfer, and the condensation steps. Eventually, a thioester linkage was formed and the protease catalytic activity was inhibited when C145 lost its ability to hydrolyze peptides.  In this conformation, the C145 residue would first be deprotonated by backbone O(H163), followed by its nucleophilic attack to the carboxylate carbonyl carbon of the clovamide molecule. The initial step of the thioesterification reaction involves protonation of the carboxylate carbonyl oxygen, which could also be assisted by possible hydrogen bonds as shown, making the group a stronger electrophile to be attacked by C145 thiolate ion. Yellow dash lines are two interatomic distances critical to the inhibitory mechanism; green dash lines represent hydrogen bonds; the magenta line indicates that the distance between C145 and the centroid of the catechol ring (5.737 Å) determines that no inhibitory reactions are favored at the catechol moiety. (B) Teo-dimensional display of non-covalent interactions between M pro and clovamide molecule docking pose 12. Numerous hydrogen bonds stabilized the conformation between the clovamide molecule and active-site amino acids, as shown in Figure 7A. Specific non-covalent interactions were laid out in the 2D diagram shown in Figure 7B. The docking study revealed a sulfur-π interaction between M49 and the catechol ring at the double bond conjugation side of the molecule at a distance, S(M49) and the centroid distance was 4.486 Å, providing extra stabilization to the protease-inhibitor complex. This sulfur-π interaction motif was reported in our previous study on inhibitory actions of Celastrol to SARS-CoV-2 main protease [45].

Clovamide and Analogs in the Zwitterionic Cys-His M pro Catalytic Site
More importantly, recent structural studies have reported the protonation states of the M pro catalytic dyad (C145 and H41) [47,48], which provided vital insights to our mechanistic studies. It was shown that, from X-ray and neutron diffraction crystallographic data, the native protease catalytic site adopted a zwitterionic (paired ionic) form, in which C145 was deprotonated in its thiolate state, hence, negatively charged, and H41 was positively charged as both of the nitrogen were protonated. Normally, the active site possesses a neutral, unreactive form, because the cysteine thiol has a pK a of 8.35 and imidazole of histidine has pK a of around 7 [49]. As a result, both species should be expected to remain non-charged under a biological pH (around 7.4). The fact that in M pro the zwitterionic form already existed prior to the substrate binding was surprising and unusual, yet critical to our analyses, since the resulting electrostatic environment had to be taken into account. In the case of the inhibitory activity in pose 12 (Figure 7), the carboxylic acid in the clovamide molecule was expected to be deprotonated under a biological pH (7.4), given that the alkyl carboxylic acid compounds had pK a values of around 4.0 to 5.0, and that the Asp and Gln free amino acids had a pK a of 4.1 [49]. Consequently, this negatively charged carboxylate group would ionically interact with the positively charged H163 residue to stabilize the protease-inhibitor complex, while providing another hydrogen bond between the H41 imidazole nitrogen and carboxylate oxygen (COO − ). This ionic interaction between the carboxylate and histidinium ion has been previously shown to stabilize the complex formed by cysteine proteases and their inhibitors; particularly, Schirmeister (2000) reported that the addition of a carboxylic acid moiety in the cysteine protease inhibitor significantly increased the rate of inhibition [50].
Before the work of Kneller et al. (2020) [48], studies suggested a ligand-free enzyme would adopt the catalytically resting state in the enzymatic activity is activated by a proton transfer from C145 to H41 imidazole upon substrate binding [51]. However, previous quantum mechanical/molecular mechanical studies have shown that this proton transfer has a high energy barrier given the measured k cat [47]. Kneller et al. (2020) [48] then reported this zwitterionic state, which allowed us to adopt this natural form in our in silico studies for the better tailoring of inhibitors.
The clovamide molecule and its three analogue species (Figure 3) was docked to this SARS-CoV-2 M pro with the zwitterionic active site (PDB ID: 7JUN). First, the docking of clovamide did not show a favorable proximity with either of the catechol rings nor other functional moieties. Docking studies of the ortho-quinone analog species showed similar enzyme-substrate conformations. Here, we described the docking of the clovamide analogue C (double-quinone clovamide) and this M pro with charged amino acids at the active site. Results of the docking double-quinone (15 poses) were analyzed and two poses, pose 3 and pose 11, were selected for the molecular dynamics cascade. Figure 8 shows the results of the dynamics cascade of pose 3. The initial docking of this pose showed C145 in close proximity to the catechol ring at the non-conjugation side of the molecule with a distance between S(C145) and the α carbon 3.821 Å; S(C145) to the centroid of this ortho-quinone ring had a distance of 4.719 Å. The dynamics cascade results showed an increase in interactions between the clovamide analogue and the catalytic site. The S(C145) to α carbon distance became 3.283 Å and its distance to centroid evolved to 4.389 Å ( Figure 8A). In addition, the native form of the C145-H41 dyad in the PDB structure 7JUN had a distance between the S(C145) and the H(H41) on the imidazole labile nitrogen of 3.007 Å. After the standard dynamics cascade, the distance became 1.985 Å, a marked decrease in distance, suggesting an electrostatic interaction between the two species, stabilizing the thiolate form of C145 assisting its binding and attack on the ortho-quinone ring. The C145 residue was expected to make a covalent linkage via a mechanism illustrated in Scheme 3(A2). The 2D diagram ( Figure 8B,C) interactions show two novel non-covalent interactions: the salt bridge and lone pair-π interactions. The salt bridge is a combination of an electrostatic attractive force (ionic bonding) and hydrogen bonding, which in this conformation happened at the Cys-His catalytic dyad. The lone pair-π interaction stabilized the substrate-enzyme complex by having the lone pair on the carbonyl oxygen of the N142 side chain interact with the π system of the orthoquinone ring. This type of interaction has been manifested in biological systems, such as its participation in the formation of DNA-protein complexes [52], and in the field of supramolecular chemistry [53]. In our case, the distance between O(N142) and the centroid of the quinone ring had a distance of 2.839 Å, stabilizing the inhibitor-enzyme complex via this lone pair-π interaction.  Figure 9A shows the outcome of the standard dynamics cascade of pose 11. In this conformation, the C145 residue was not close to either of the quinone rings but was in proximity to the amide moiety towards the center of the molecule. As shown in Figure 9A, the S(C145) had a short distance with the carbonyl carbon of the amide group, 3.158 Å. The catalytic C145 residue was stabilized by non-covalent interactions provided by (1) attractive ionic forces with positively a charged H41 residue, with a H(H41) to S(C145) distance of 2.032 Å after the dynamics study. The initial docking kept the native distance of 3.007 Å as shown previously; (2) hydrogen bonding with the carboxylic acid moiety with a H(COOH) to S(C145) distance of 2.222 Å. The docking pose did not manifest this interaction with a long distance between H(COOH) and S(C145) of 6.587 Å. With these stabilizing forces, we proposed that the C145 would attack the carbonyl carbon of the amide group and perform inhibition by the previously described mechanism that serine/cysteine proteases adopt when they hydrolyze peptide chains [54]. Upon substrate binding, the thiolate ion attacked the carbonyl carbon forming a negatively charged oxygen anion species in a tetrahedral intermediate, which was stabilized by hydrogen bonding with the polypeptide backbone N-H, also known as the oxyanion hole in the serine/cysteine protease catalytic mechanism. In our case, similarly, this oxyanion hole would be provided by G143 and S144 residues as shown in other studies as well. This oxyanion hole region was labeled in Figure 9A as the dotted orange circle. Green dashes between carbonyl oxygen with backbone N-H of G143 and S144 highlighted this ionic stabilization. The following steps were proposed to be similar as found in the Chymotrypsin cleavage mechanism, in which acyl enzyme intermediates were formed followed by the leaving of the first amine product. Next, the approximation of a water molecule led to the formation of the second tetrahedral intermediate, which resulted in the formation and departure of a second acid product.
We proposed that, by acting to cleave the peptide bond of the amide moiety in the clovamide molecule, the catalytic activity of the C145 residue could be inhibited. However, unlike the canonical serine/cysteine protease mechanism, the clovamide in our case was not the substrate; rather, it acted as the inhibitor and we expected it to be an irreversible inhibitor. A difficulty associated with this hypothesis was the fact that upon the completion of the peptide cleavage, the active site was eventually reformed as the catalytic cycle completed. This would not serve our purpose as we anticipated a covalent inhibition mechanism. Despite this fact, we argue that this mechanism actually helped the formation of a neutral Cys-His dyad at the end of the acting cycle, creating barriers for the catalysis to be activated. The final step of this mechanism was expected to be the formation of a neutral C145-H41 catalytic dyad with a protonated C145 (S-H) and deprotonated H41 both remaining uncharged. To activate the catalytic activity, the initial step had to involve the transfer of this proton to its zwitterionic state, an essential first-step with a great energy barrier as explained previously [48]. We, therefore, argue that this mechanism inhibited the catalytic activity not by covalent modification (although covalent bonds were formed at intermediate steps), but primarily by creating a barrier for catalysis by forming the neutral catalytic dyad. These specific steps of series of reactions are illustrated by Scheme 5.
Other non-covalent interactions that stabilized the inhibitor-enzyme complex are also shown in Figure 9B,C. A novel non-covalent interaction shown here was the π-anion interaction. This interaction was similar to the earlier introduced lone pair-π interaction and played a key role in the stability of the overall complex [55]. Figure S4 provides a zoomed-in picture of this interaction.
The further examination of the reactivity of ortho-quinone clovamide species was performed by analyzing their molecular orbitals. An energy optimization of clovamide double-quinone species (derivative C, Figure 3) showed LUMO orbitals concentrated at the ortho-quinone ring moiety, suggesting reactivity at the quinone ring ( Figure S5). Specifically, a preferred reactivity at the quinone near the double-bond conjugation side was suggested, which was the case both for the natural clovamide molecule ( Figure S5A) and its anionic form as the carboxylic acid group was deprotonated ( Figure S5B). This preference was in accordance with our docking studies performed on quinone clovamide derivative C, as shown in Figure 6. The proton transfer step with high energetic barriers which prevents the catalytic site to reform its charged, zwitterionic state that is essential for its peptide hydrolysis catalytic activity.

Allosteric Inhibitions: Clovamide Acts as SARS-CoV-2 M pro Non-Competitive/Mixed-Type Inhibitor
Last, we explored potential inhibitory activities by clovamide on the allosteric sites of SARS-CoV-2 M pro . Most research projects related to COVID-19 drug design targeting the main protease have been heavily focused on the catalytic domain. Recent studies [56,57] explored small-molecule drugs that bind to the main protease allosteric sites and suggest that they act as non-competitive inhibitors. Specifically, Gunther et al. (2021) [57] performed high-throughput X-ray Crystallographic studies and identified multiple organic small molecules that bind to two distinct allosteric sites, respectively, at domain II and domain III of the protease monomer, with strong antiviral activities.
We, respectively, performed molecular docking studies on the two allosteric sites as reported [57]. Figure 10 shows the positions of the two allosteric binding sites in the protease dimer: site one was located at a three helices bundle region of domain III, in close contact with the opposite protomer (binding pocket shown in red, Figure 10B); site two was located in between domain II and domain III, nonboring to site one (binding pocket shown in yellow, Figure 10A). The inhibitors with potent antiviral activities that were found in the protease-inhibitor crystal structures were spatially overlapped with the clovamide conformations from our docking results. At both allosteric sites, the overlaps between clovamide and known strong inhibitor molecules were promising. Allosteric site one located at the interface of the two homodimers and, therefore, any conformational changes as a result of inhibitor bindings would lead to changes in the opposite dimer. We set docking studies at this site (apo form of the PDB 7AXM structure) and within the returned 20 poses, pose 8 showed an interesting spatial arrangement that allowed clovamide to interact with site one residues through every part of the molecule ( Figure 11A). First, the top portion of the clovamide molecule ( Figure 11B) bound into the hydrophobic pocket, primarily due to hydrophobic interactions by the aromatic catechol ring and non-polar amino acid residues (I213, L253, Q256, V296, and V297). The two hydroxy groups on the catechol ring were pointed outwards, facing the solvent-exposed side and away from the pocket. We proposed that, similar to what was found in the complex structure with inhibitor Pelitinib [57], the binding of clovamide would cause a conformational change of the Q256 residue, swinging away from the hydrophobic pocket ( Figure 11B). Lastly, the carboxylic acid moiety also pointed outward and faced away from the allosteric site, which was largely hydrophobic. This binding pose had a binding energy of −5.1 Kcal/mol. Here, we used the parameter of binding energy (also known as binding affinity) as part of the docking outcomes to indicate the strength of interactions between the enzyme and the inhibitor: the more negative the energy, the stronger the interactions; hence, the more stable the reversibly bound complex was. The lower portion of the molecule was in close contact with the opposite protomer, in a region that was close to the catalytic active site. Figure 11C depicts the binding mode; the lower part of clovamide was oriented towards the substrate binding pocket. Here, a short peptide was included (extracted from the crystal structure of the a protease-substrate complex of the highly homologous SARS-CoV M pro , PDB ID: 2Q6G) [58] to indicate the substrate binding sites around the catalytic domain. The aromatic ring of clovamide (shown in cyan) was in proximity to the Q6 residue, which was the first residue after hydrolysis P1 by the catalytic dyad (TSAVLQ↓SGFRK, ↓ indicates the hydrolysis site). This was a critical pocket (S1 pocket) that directly involved in the catalytic reactions at the active site whose integrity was indispensable for enzymatic activity [59]. In addition, similar to that found in the complex crystal structure of Pelitinib, we argue that clovamide binding could also cause a conformational change of the N142 residue ( Figure 11C) and, therefore, affect the integrity of the active site and substrate binding pockets as well. As a result, both the catalytic domain active site and the substrate binding sites could be collapsed and the catalytic activity of the protease be inhibited.
Allosteric site two was located in between domain II and III of the monomer form, and was not directly in contact with the opposing protomer. However, the integrity of this allosteric site was crucial in ensuring normal catalytic functions of the protease. Our docking studies suggested a great binding affinity towards this site: the best pose had a binding energy of −8.3 Kcal/mol and is shown in Figure 12A. Similar to the allosteric site one, the upper portion of the molecule buried into a hydrophobic pocket consisted of residues Q107, P108, G109, Q110, V202, I249, and T292. The peptide (amide) moiety at the center of the molecule interacted with its nonboring aromatic F294 residue via an amide-π interaction [60]. A critical interaction relied on the carboxylic acid group of clovamide molecule; the carboxylic acid moiety established strong hydrogen bonding with the D153 residue as well as the R298 residue. As a result of the COOH-D153 hydrogen bond, the D153 residue underwent spatial rearrangement and swung towards the COOH group. The C-α of D153 migrated 1.490 Å towards the inhibitor [57]. This conformational change led to the subsequent movement of the nonboring Y154 and D155 residues ( Figure 12B); in particular, Y154 underwent a dramatic rearrangement that caused the loop 153-155 to migrate inwards, with a C-α movement of 2.819 Å [57].
Most importantly, the carboxylic acid group was able to form hydrogen bond and ionic interactions with the side chain of R298 residue, a residue that was shown to be essential for the protease dimerization process. R298 acted as a "switch" that directly controlled the conformations of the N and C termini of the protein; any changes or mutations would lead to a switch from a dimer to monomer, completely eliminating the catalytic activity [61], since the M pro is only catalytically active when it is a dimer [6]. In reference [57], the inhibitor AT7519 showed antiviral activity and was bound to this allosteric site [57]. We, therefore, argue that clovamide also had a favored moiety to allow similar inhibitor-residue interactions and, hence, inhibit the protease catalysis by hindering the dimerization when R298 is engaged in salt-bridge interactions.
Furthermore, we performed molecular docking studies to explore the potentials of other functional groups in engaging interactions of R298 residue. The carboxylic acid group was substituted by carboxylate, a hydrogen atom, a hydroxy group (Figure 12E), and a methyl group. When it was substituted with a carboxylate (COO -) group ( Figure 12C), the strong hydrogen bond remained, with a O(COO -) to N(R298) heteroatomic distance of 2.808 Å. Replacing the carboxylic acid group with a hydrogen atom, as expected, eliminated the interaction; however, we observed a reorientation of the lower catechol group from its original position (close to N151 as shown in Figure 12A) towards the R298 residue, forming a hydrogen bond with a heteroatomic distance of 3.109 Å ( Figure 12D). The hydrogen bonding interactions between R298 and the chiral group could be reformed by introducing groups available to form hydrogen bonds, such as a hydroxy group (Figure 12E), despite a weaker bond (3.227 Å). Interestingly, no favored interactions could be observed when the group was substituted with a methyl group. We argue that this may be due to steric effects given that a methyl group is bulkier than hydrogen, creating a steric hindrance even when the catechol moiety attempts to reorient and forms hydrogen bonds such as the hydrogen-substituted species. Therefore, future rational drug design should take into account steric factors since, in the case of clovamide, its aromatic group on the top favored the hydrophobic pockets; thus, limiting the structural freedom of the molecular binding conformations.

Discussion
As of today, no antiviral drugs have been officially approved by the World Health Organization (WHO) for COVID-19 therapeutics [62]. Prior candidates of small-molecule pharmaceuticals, such as hydroxychloroquine and remdesivir [63], have been shown clinically to be unsuccessful in the prevention of COVID-19. In fact, the use of these two drugs was recommended against by the WHO [64]. A recent study showed that underlying the supposed antiviral activities of hydroxychloroquine were actually a result of a shared mechanism of phospholipidosis, which depends on physiochemical properties rather than expected target-based antiviral effects [65]. It is, therefore, of interest to further the efforts of developing new antiviral drugs to assist in combating the pandemic. The SARS-CoV-2 M pro is a critical enzyme in the viral replication cycle; hence, a promising target for antiviral agents. In searching for drug lead molecules, natural products provide a scaffold and insights as drug precursors. Extensive studies suggest that natural compounds, especially phenolic compounds, are capable of inhibiting SARS-CoV-2 M pro ; hence, they may act as COVID-19 pharmaceuticals [13][14][15][16][17]. This work focused on a natural antioxidant, clovamide, found in the phenolic fraction of cacao, whose extracts have been shown to be able to inhibit the M pro activity [18]. Our computational studies suggested mechanisms by which clovamide and its derivatives inhibit M pro when bound to catalytic and allosteric sites of the enzyme. These studies also revealed important structural skeletons for designing M pro inhibitors.
Initial docking studies of clovamide and apo protein (PDB ID: 6LU7) did not show favored conformations that suggested inhibition via the proposed Michael addition mechanisms (Scheme 1). No reactivity was shown with the α, βunsaturated amide moiety. However, our subsequent docking and dynamics studies of clovamide quinone derivatives and the protease gave us promising results in which the C145 catalytic residue was in proximity to the quinone centroid (3.450 Å), indicating chemical reactivity. Three conjugate addition mechanisms were proposed for the irreversible, covalent inhibition of clovamide quinone derivatives on the catalytic site of M pro . The mechanisms (Scheme 3) were rendered likely by experimental evidence: upon radical scavenging (anti-or pro-oxidant effects), clovamide becomes oxidized to ortho-quinone species [43], which, then, readily react with thiol groups to form thiol adducts [42], forming the linkage between the sulfur atom of C145 and carbon on the quinone ring. More importantly, these quinone derivatives of clovamide can be naturally found in cacao as a result of both non-enzymatic and enzymatic oxidation mediated by PPO, which also involves a subsequent crosslinking with proteins [26,27]. In addition, clovamide has been determined to be a broad-spectrum enzyme inhibitor, shown to inhibit proteases by a generic mechanism involving oxidation to quinone and the formation of covalent linkages with reducing groups such as thiol [26]. These findings greatly support our mechanistic proposal and complement our computational results.
Clovamide is a versatile natural compound with various functional groups to exploit in enzymatic inhibition processes. Thus, we also performed computational studies to examine the role of other functional groups in inhibiting the catalytic activity of the protease. The carboxylic acid group is a chiral moiety in the clovamide molecule, able to provide binding affinity via hydrogen bonding. Our exploration also suggested the possibility of carboxylic acid acting as a chemical warhead that provides chemical reactivity, such that the deprotonated thiolate C145 could attack the carbonyl carbon. We proposed this condensation mechanism (Scheme 4); yet, there lacked experimental support for this reaction mechanism. Moreover, a change in the binding mode was observed when we docked clovamide against the active site of the protease with a zwitterionic catalytic dyad (PDB ID: 7JUN) [48]. The docking model suggested reactivity between the native deprotonated thiolate group and the carbonyl carbon of the α,βunsaturated amide moiety. Considering the similarity of this conformation with the initial conformation when cysteine/serine proteases catalyze peptide substrates, we proposed a similar mechanism in which a clovamide molecule binds and becomes cleaved. We argue that this reversible binding causes the catalytic dyad to reform its neutral state from the native zwitterionic state. As a result, the catalytic site must overcome a high energy barrier to achieve the deprotonation of C145 to start catalysis. This mode of inhibition, however, is unprecedented and experimental data are required to validate the proposal.
The allosteric inhibition of SARS-CoV-2 M pro can also be achieved by clovamide. We docked the clovamide molecule to the two allosteric sites of M pro discovered by highthroughput X-ray screening studies [57]. We selected the enzyme complexes with inhibitors showing the strongest antiviral effects and docked clovamide against the apo structures of those known complexes. At both sites, clovamide showed great binding affinities with strong non-covalent interactions. At site one, clovamide inhibited protease activity by extending itself to the substrate binding pocket of opposing protomer and affecting the integrity of that active site; at site two, clovamide bound tightly and formed saltbridge interactions with the R298 residue which is the switch that controls the protease dimerization processes.
The availability of high-resolution crystal structures of SARS-CoV-2 M pro and in complex with known inhibitors enabled us to use computational approaches for preliminary structure-based studies. However, clovamide will need to be obtained and tested against its antiviral efficacy. Synthetic routes of clovamide and a variety of its analogues with catechol substituted with different R groups have been reported [66,67]. The synthetic schemes involve as simple as a two-step synthesis to obtain clovamide products [66].
Last, we propose the following future experimental assays to test our computational models of clovamide as a COVID-19 drug candidate. First, performing kinetic assays to assess the efficacy of clovamide acting as a SARS-CoV-2 main protease inhibitor. Kinetic information is essential for evaluating enzymatic inhibition. By obtaining the rate constant data and establishing Michaelis-Menten double-reciprocal plots, one could deduce what type of inhibitor clovamide is and gain insight regarding possible steps in the inhibition process [68]. Secondly, obtaining structural data to elucidate the inhibitory mechanism. A crystal structure of the protease-inhibitor complex would directly confirm which of our proposed inhibitory mechanism(s), either catalytic or allosteric, is favored. Once a structure is solved, a detailed analysis of intermolecular interactions between the inhibitor and the protease could also suggest critical features of the inhibitor that is important for binding and reactivity (structure-activity relationships). Finally, conducting cell-based assays to examine the antiviral potential of clovamide. For instance, PCR techniques could be used to assess antiviral effects of clovamide when applied to SARS-CoV-2-infected Vero cells [6].

Conclusions
Our work suggested that clovamide may be a promising COVID-19 drug lead molecule that targets the SARS-CoV-2 main protease. We hereby show that clovamide is a versatile natural compound with various functional groups, providing both binding affinities and chemical reactivities to inhibit SARS-CoV-2 main protease activities. Our molecular docking and dynamics studies revealed possible reaction mechanisms by which the oxidized ortho-quinone moiety readily reacts with C145 thiolate residue as an irreversible inhibitor. We propose that this ortho-quinone moiety is a favored group providing warhead reactivity in inhibiting cysteine proteases. When the catalytic dyad adopts a native zwitterionic state, we propose that clovamide could bind and prompt the neutral state to reform, creating a high-energy proton transfer barrier. We also suggest that, by reversibly binding to the two allosteric sites of the protease, clovamide is also capable of inhibiting the enzyme by inducing local enzymatic conformational changes and altering bonding environments of critical residues; hence, inhibiting enzymatic activity by allosteric inhibition mechanisms. We hope this work promotes future experimental studies to test our computational models. Eventually, insights from this work will benefit a future structurebased rational design of SARS-CoV-2 main protease inhibitors and future developments of COVID-19 pharmaceuticals.