In silico approach identified benzoylguanidines as SARS-CoV-2 main protease (Mpro) potential inhibitors

Abstract The coronavirus disease-2019 (COVID-19) pandemic, caused by the novel coronavirus severe acute respiratory syndrome-coronavirus-2 (SARS-CoV-2), became the highest public health crisis nowadays. Although the use of approved vaccines for emergency immunization and the reuse of FDA-approved drugs remains at the forefront, the search for new, more selective, and potent drug candidates from synthetic compounds is also a viable alternative to combat this viral disease. In this context, the present study employed a computational virtual screening approach based on molecular docking and molecular dynamics (MD) simulation to identify possible inhibitors for SARS-CoV-2 Mpro (main protease), an important molecular target required for the maturation of the various polyproteins involved in viral replication. The virtual screening approach selected four potential inhibitors against SARS-CoV-2 Mpro. In addition, MD simulation studies revealed changes in the positions of the ligands during the simulations compared to the complex obtained in the molecular docking studies, showing the benzoylguanidines LMed-110 and LMed-136 have a higher affinity for the active site compared to the other structures that tended to leave the active site. Besides, there was a better understanding of the formation and stability of the existing H-bonds in the formed complexes and the energetic contributions to the stability of the target-ligand molecular complexes. Finally, the in silico prediction of the ADME profile suggested that LMed-136 has drug-like characteristics and good pharmacokinetic properties. Therefore, from the present study, it can be suggested that these structures can inhibit SARS-CoV-2 Mpro. Nevertheless, further studies are needed in vitro assays to investigate the antiviral properties of these structures against SARS-CoV-2. Communicated by Ramaswamy H. Sarma


Introduction
Coronaviruses (CoVs) belong to subfamily Coronavirinae, from family Coronaviridae, first noticed in 1947 and divided into four genera: a, b, c, and d. A novel coronavirus of the b genus, known as severe acute respiratory syndrome-coronavirus-2 (SARS-CoV-2), emerged in China in December 2019 and later was identified as the coronavirus disease-2019 (COVID-19) etiologic agent. A few months later, this disease was declared a pandemic in March 2020 by World Health Organization (WHO) (Hage-Melim et al., 2020). People infected not only can develop mainly symptoms such as fever, cough, and difficulty in breathing, but also could have an asymptomatic aspect. In this way, SARS-CoV-2 is reported to cause respiratory problems, besides being recently described as also causing neurological problems (Hagar et al., 2020). Until, the present moment (August 31, 2022), this virus has been responsible for more than 599.825.400 confirmed cases of infected people and 6.469.458 deaths worldwide (World Health Organization, 2022a) The scientific and medical community worldwide has been working hard to find a solution to control the pandemic, urgently developing and approving vaccines that accounted for 12.449.443.718 people immunized (August 31, 2022), (World Health Organization, 2022a). The rapid development of vaccines has been crucial in combating the COVID-19 pandemic. However, some issues are complex to immunizing the global population, such as anti-vaccination campaigns and contraindications for people with individual health problems (Owen et al., 2021). One of the biggest challenges of vaccination is global distribution and access to the immunizer, given that the most significant part of vaccine doses is concentrated in a small group of rich countries, which administer about three-quarters of all vaccine doses. Until August 2021, six of the ten countries with the highest number of deaths per capita from COVID-19 (Tunisia, Georgia, Botswana, Eswatini, Namibia, and South Africa) had less than 10% of their populations fully vaccinated. The other four (Fiji, Malaysia, Cuba, and Kazakhstan) have fully vaccinated more than a third of their population (Schaefer et al., 2021). Consequently, vaccine inequality can lead to possible new endemic foci and the emergence of new variants. Therefore, the need for an easily accessible antiviral drug is evident. Faced with the need for therapies to treat COVID-19 in urgent care, the US Food and Drug Administration (FDA) approved, in October 2020, the emergency use of the antiviral drug Remdesivir (Beigel et al., 2020). Remdesivir is a phosphoramidate prodrug that is metabolized in cells to produce the active form that acts as a nucleoside analog, being incorporated by the RdRp into the growing RNA product, inhibiting RNA-dependent RNA polymerase (RdRp) impairing the viral replication of coronaviruses, including SARS -CoV-2 (Kokic et al., 2021). Although this drug has been recommended for use in adult and pediatric patients 12 years of age with a minimum weight of 40 kg who require hospitalization, it still requires substantial evidence of efficacy and safety for the intended purpose, in addition to the high cost of implementing of the therapy (FDA, 2021a).
Furthermore, the FDA authorized in December 2021 the emergency use of Paxlovid produced by Pfizer (FDA, 2021b). Paxlovid is a combination of nirmatrelvir and ritonavir. The first drug inhibits viral replication by inhibiting the major protease (M pro ) by forming a reversible covalent thioimidate adduct between the nitrile group of nirmatrelvir with the catalytic residue of Cys145, and the second one is responsible for delaying the degradation of nirmatrelvir (Owen et al., 2021). Therefore, treatment is recommended for adult and pediatric patients (12 years or older, weighing 40 kg or a minimum of 88 pounds) at a high risk of progression to severe COVID-19, including hospitalization or death (FDA, 2021b). In addition, in April 2022, the Guideline Development Group (GDG) published a new recommendation for Paxlovid (nirmatrelvir and Ritonavir), indicating it for patients with non-severe COVID-19 who are at increased risk of hospitalization (World Health Organization, 2022b).
In addition to Remdesevir and Paxlovid, other therapeutic options are available, such as the antiviral molnupiravir; IL-6 receptor blockers like tocilizumab or sarilumab in combination with corticosteroids; Janus kinase inhibitors (JAK), specifically baricitinib, an oral drug that suppresses overstimulation of the immune system, recommended for patients with severe or critical COVID-19, administered with corticosteroids. In the absence of JAK inhibitors (baricitinib) and IL-6 receptor blockers (tocilizumab or sarilumab), doctors may prescribe two other JAK inhibitors (ruxolitinib or tofacitinib). In addition, the monoclonal antibody sotrovimab is recommended to treat mild or moderate COVID-19 in patients at high risk of hospitalization. Finally, casirivimab-imdevimab, intended for patients with seronegative status, where rapid viral genotyping is available and confirms infection with a susceptible SARS-CoV-2 variant, represents an alternative to sotrovimab (World Health Organization, 2022a). Despite containing a wide range of drugs for the treatment of COVID-19, few directly target the SARS-CoV-2 virus, demonstrating that as new variants emerge, the development of direct-acting antivirals for outpatient treatment is necessary for COVID-19.
The drug development process is a significant challenge in the pharmaceutical industry. Drug discovery includes target selection, hit identification, lead optimization, and preclinical and clinical studies that require substantial time and money (Maia et al., 2020). In order to minimize the cost and time for these processes, computational tools are applied in researching new bioactive compounds. A common strategy is the virtual screening technique of extensive chemical libraries by molecular docking into attractive molecular targets (Kontoyianni, 2017). The non-structural and structural proteins are promising targets for designing and developing antiviral agents against SARS-CoV-2. In this context, the nonstructural main protease (M pro ) is conserved among coronaviruses and is responsible for cleaving polyproteins at 11 different sites to produce shorter non-structural proteins crucial for viral replication (Jin et al., 2020). Thus, inhibitors interacting with this target can be great therapeutic or antiviral agents against COVID-19 (Ghosh et al., 2021). Therefore, given the evident need for new antiviral therapies and considering the importance of M pro in viral replication, several studies have evaluated different sets of molecules to find lead compounds (Arif, 2022;Lokhande et al., 2021;Sharma et al., 2022;Yu et al., 2021). In this context, the present work was carried out to identify possible new inhibitors of SARS-CoV-2 M pro from our in-house library (LaSMMed Chemical Library-LMCL), which contains a diversity of 313 structures from different chemical classes, such as benzoylguanidines, benzoylthioureas, benzoylureas, carvacrol derivatives, coumarins, hydantoins, indoles, thiazolidinones, thiohydantoins, and oxazolidinones. These scaffolds can be considered privileged structures (DeSimone et al., 2004) and present several in vitro biological activities, mainly as antimicrobial and antiprotozoal agents (Benatti et al., 2020;Brito et al., 2020;Camargo et al., 2022;de Carvalho et al., 2019;Pereira et al., 2021;Santiago-Silva et al., 2021). In addition, there are reports of substances with these scaffolds that showed good antiviral activities against different RNA viruses (Ghosh et al., 2018;Hishiki et al., 2017;Hwu et al., 2019;Pilau et al., 2011;Tijsma et al., 2016;Xu et al., 2020). Therefore, given previous evidence of antiviral activity, we decided to exploit our in-house library (LMCL) as chemical space in a simplified and rapid in silico protocol to identify promising hits for discovering new treatments against COVID-19 (Fig. 1).

Protein structure and LMCL
The study was performed with crystallographic structures of two M pro , one of them in complex with an N3 inhibitor (PDB ID: 6LU7, chain A, resolution 2.16 Å) (Jin et al., 2020) and the other with an inhibitor a-ketoamide (13 b) (PDB ID: 6Y2F, chain A, resolution 1.95 Å) .
The missing residues of protein structures were added using the CHARMM-GUI platform (http://www.charmm-gui. org/) (Lee et al., 2016). After, water molecules and co-crystallized inhibitors were removed from the structure.
The 3D structures from the LMCL were drawn using ChemDraw and converted to pdbqt using the Open Babel (O'Boyle et al., 2011). Finally, geometry optimization was performed using the MM2 force field implemented by ChemBio3D v.12.0 program (PerkinElmer Informatics) (Evans, 2014).

Molecular docking-based virtual screening
Virtual screening experiments were performed using AutoDock Vina 1.1.2 (Trott & Olson, 2009)   was centered on co-crystallized inhibitors: M pro ID: 6LU7 (x: À 10.729; y: 12.417; and z: 68.816) and for the M pro ID: 6Y2F (x: 10.879; y: À 0.251; and z: 20.753). The box size for both (25 � 25 � 25 Å 3 ) was determined to cover the active site. The grid box to Autodock4.2 was centered in the same coordinates previously described to Autodock Vina, but with a grid dimension of 60 � 60 � 60, and the spacing was 0.375 Å. The Lamarckian Genetic Algorithm (LGA) was used to perform the simulation. Both software returned 10 binding poses per compound, and the highest-scoring poses were used to calculate the arithmetic mean of the affinity energy in both enzymes. The root-mean-square deviations (RMSD) calculations were carried out using PyMOL software (DeLano, 2002).

Molecular dynamics simulations
Molecular dynamics (MD) simulations were carried out using GROMACS 2018.1 package (Abraham et al., 2015;Berendsen et al., 1995) with CHARMM36 force field (Huang & MacKerell, 2013) from the best pose obtained by virtual screening. In addition, the enzyme topology preparation was carried out according to described by our group (Albuquerque et al., 2020;Lima et al., 2015).
The systems were submitted to energy minimization steps using the steepest descent method with a convergence criterion of 1000 kJ mol À 1 nm À 1 , followed by gradient conjugate until the power 100 kJ mol À 1 nm À 1 .
The minimized systems were submitted to two consecutive equilibration steps, considering 300 K of temperature and 1 bar of pressure. Both steps were realized during 1 ns with position restraint for enzyme and ligand. The first step considered the number of particles, volume, and temperature constant (NVT ensemble), and the second step considered the system as isothermal-isobaric (NPT ensemble).
The temperature control was realized with a V-rescale thermostat (Bussi et al., 2007) and the pressure control with a Parrinello-Rahman barostat (Parrinello & Rahman, 1981). First, the hydrogen atoms in the system were frozen using the LINCS algorithm (Hess et al., 1997;Linear Constraint Solver). Next, the long-distance electrostatic interactions were treated using the PME algorithm (Particle-Mesh Ewald) (Darden et al., 1993;Essmann et al., 1995), and the radius cut-off was applied to the van der Waals, and Coulomb interactions were 1 nm. After the equilibrium, MD simulations were performed considering the NPT set without any position restriction, using two fs of integration time for 150 ns.
All M pro -ligand complexes were evaluated in RMSD and root-mean-square fluctuation (RMSF) using gmx rms and gmx rmsf modules available in the GROMACS package. In addition, the hydrogen bonds (H-bond) analysis was evaluated with the gmx hbond module considering the cut-off distance H-bonds (D-H … A) until 0.40 nm and the cut-off H-D-A angle until 30 � .
The clustering analysis was carried out with the gmx cluster module using the gromos method, an RMSD cutoff of 0.2 nm from the M pro -ligand complex. After that, a total of 200 frames of the most representative cluster were used to calculate the binding free energy (DGbind) (Albuquerque et al., 2020;. The DGbind to the M pro -ligand complex was calculated using the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method through the g_MMPBSA module version 5.1.2 (Kumari et al., 2014).
The simulation trajectories were visualized with VMD v.1.9.3 software (Humphrey et al., 1996), and the RMSD time series were plotted as a function of time with the XMgrace tool (Turner, 2005). The H-bond frequency was calculated with the software HbMap2Grace using a cut-off radius of 4.0 Å (Gomes et al., 2009).

Molecular docking-based virtual screening
Firstly, we performed a virtual screening against the M pro from SARS-CoV-2, a conserved coronaviruses enzyme essential for viral replication (B� aez-Santos et al., 2015), to investigate the ability of the structures to fit into the protein binding pocket and interact with Cys-His catalytic dyad (His41 and Cys145) (Jin et al., 2020). Compounds that specifically bind His41 and Cys145 can weaken or block the catalytic activity of M pro (Mengist et al., 2021). Proteins are flexible and can be crystallized into different conformations, resulting in side chains variations at the active site (Stoddard et al., 2020). Therefore, we decided to use two crystal structures of SARS-CoV-2 M pro ID: 6LU7 (Jin et al., 2020) and the M pro ID: 6Y2F (Zhang et al., 2020), which have side chains at the active site in different conformations, to evaluate as the diversity of accessible variations can produce different binding poses for a possible inhibitor.
We performed the docking protocol validation through the redocking simulation method to the ligands N3 and 13 b co-crystallized with 6LU7 and 6Y2F complexes, respectively. The pose results with the highest scoring for both enzyme structures show RMSD values of 1.41 and 1.87 Å for the 6LU7 and 6Y2F (Figure 2), respectively, which are in an acceptable range (Hevener et al., 2009).
After the protocol validation of molecular docking, we accomplished a virtual screening to investigate the LMCL, which contains a structural diversity distributed in ten main chemical classes: benzoylguanidines, benzoylthioureas, benzoylureas, carvacrol derivatives, coumarins, hydantoins, indoles, thiazolidinones, thiohydantoins, and oxazolidinones ( Figure 3). Among these classes, benzoylguanidines could have three possible amino-imino tautomeric forms (O'Donovan et al., 2013), and the thiazolidinones were obtained as racemic mixtures. Therefore, we decided to  evaluate the influence of tautomerism and stereochemistry on the affinity with M pro , yielding a total of 313 structures. These privileged structures were designed using medicinal chemistry strategies, such as molecular hybridization, conformational restriction, and bioisosterism. Moreover, they presented several in vitro biological activities, mainly antimicrobial and antiprotozoal agents (Benatti et al., 2020;Brito et al., 2020;Camargo et al., 2022;de Carvalho et al., 2019;Pereira et al., 2021;Santiago-Silva et al., 2021). Besides, in the literature, many reports about the antiviral potential of these scaffolds, mainly against RNA viruses, were found (Cihan- Ust€ unda� g et al., 2016;Harnden et al., 1978;Musella et al., 2016;Pehlivan et al., 2018;Sanna et al., 2018).

€
Thus, after the virtual screening simulations, the 313 structures showed average binding affinities ranging from À 8.25 to À 3.9 kcal/mol. We then selected four structures as possible inhibitors of SARS-CoV-2 M pro since their predicted average binding affinities were below À 8 kcal/mol (Bepari & Reza, 2021). Based on our analysis, the four selected structures are benzoylguanidines, LMed-136, LMed-137, an LMed-136 tautomer, LMed-110, and LMed-117. (Table 1). Afterward, we performed molecular docking for the best four structures from our virtual screening. These analyses confirmed that these four structures could occupy the active site of the main protease of SARS-CoV-2. In addition, the four structures showed similar binding modes despite using two proteins (6LU7 and 6Y2F) with side chains having different conformations at the active site (Figures 4 and 5).
Analyzing the 6LU7-LMed-136 (Figure 4(a)) complex was stabilized by a more significant number of H-bonds interactions, six in total. In addition to the interactions involving the guanidinium group and Glu166 (S4), Gln189 (S4), and His 164 (S1) previously observed, also three other interactions are present, such as the nitro group and Thr190 and Gln192, both at S4 subsite and, finally, the carbonyl's benzoyl group with the catalytic His41 (S1'). However, the 6Y2F-LMed-136 complex ( Figure 5(a)) showed H-bonding interactions only with the amino acid residues His164 (S1), Thr190 (S4), and the catalytic histidine (His41) at the S1' subsite. Inversely, to the 6LU7-LMed-137 complex (Figure 4(b)), there was a decrease in the H-bonding interactions compared to the complexes of its tautomer LMed-136 (Figure 4(a)), being four and six H-bonds, respectively. Furthermore, these Hbonds occurred only between the nitro group with Tyr54 (S2) and His164 (S1) and the guanidinium group with Glu166 and Gln189 at the S4 subsite. In comparison, the main alteration in the 6Y2F-LMed-137 complex ( Figure 5(b)) was the more one loss of the H-bond, this time, with Gln189 (S4).
For the 6LU7-LMed-110 protein-ligand complex ( Figure  4(c)), we observed four H-bonds: two involving the guanidinium group with the Gln189 side chain at the S4 subsite; one interaction of this same group with the carbonyl of the backbone of His164 at the S1 subsite, and the last one between the carbonyl of benzoyl group with the nitrogen of the backbone of Glu166 at the subsite S4. In contrast, the main change observed in the 6Y2F-LMed-110 complex (Fig.  5C) was the loss of the two H-bonding interactions with Gln189 residue (S4). For LMed-117, we observed H-bonds with the similar interactions to the LMed-110, that is, with the amino acid residues His164 (S1), Glu166 (S4), and Gln189 (S4) in the 6LU7-LMed-117 (Figure 4(d)) complex, and only with residues His164 (S1), and Glu166 (S4) in the 6Y2F-LMed-117 ( Figure 5(d)) complex.
In summary, the docking analyses indicated that the four selected structures have a greater affinity for the S4 subsite, performing essential hydrogen bonds that can influence the enzymatic activity of M pro . Furthermore, in their study, Fu et al. observed that Boceprevir exhibits a mode of interaction similar to our study with the S4 subsite (Fu et al., 2020).

Molecular Dynamics simulations
We performed the molecular dynamics simulation (MDS) during 150 ns to evaluate the behavior of the four promising compounds (LMed-136, LMed-137, LMed-110, and LMed-117) in an aqueous system. Firstly, we evaluate the RMSD profiles for a-carbon atoms (RMSD-Ca) of M pro . The results suggested a significant shift compared to the docking complex. However, the protein structure stabilizes after the first nanoseconds of simulation with a standard deviation below 2 Å (Table 2 and Figure S1, supplementary material).
RMSD analysis of inhibitors bound to M pro reveals an orientation change at the binding site compared to the docking pose. With the LMed-136 (3.53 ± 1.49 Å) and LMed-137 (8.29 ± 3.06 Å), it is possible to observe sudden changes in orientation at the binding site compared to the anchoring pose or indicating a tendency to leave the active site ( Albuquerque et al., 2020;de Souza et al., 2019;Uelisson da Silva et al., 2022 ). About LMed-110 (1.86 ± 0.37 Å) moves below 2.0 Å, demonstrating a slight change in orientation, but returns to the starting position. Despite the LMed-117 (4.43 ± 0.40) moving above 2.0 Å, after 5 ns, there is a tendency for the pose to stabilize at the active site (Fig. 6).
Besides, the mobility of residues was evaluated through the RMSF analysis. To facilitate the analysis of residue mobility during MDS, we identified which residues are close to inhibitors up to 5 Å. Therefore, this analysis (Fig. S2) indicates that the significant fluctuations (> 3 Å) are more than 5 Å away from the active site, suggesting that the presence of these inhibitors does not affect the stability of the enzyme (Albuquerque et al., 2020).
The analysis of the H-bonds of the M pro complexes with LMed-136,  was performed during the 150 ns of MDS, since it is a crucial parameter to predict the binding mode of the ligands in the active site of M pro , as these interactions play an important role in directing the binding of an inhibitor to the receptor and also contribute in an important way to the binding affinity (Yunta, 2017). Analyzing the complex with the ligand LMed-136 and its tautomer LMed-137, we see that both perform more H-bond interactions than the ligands discussed above, but the bonds have a short lifetime. For example, the ligand LMed-136 ( Fig 7A) performed two H-bonds with Gln189, one with 42.02% and the other with a short lifetime (1.58%). The other short lifetime H-bonds are two with the residue Gly143 (1.02% and 1.06%) and two with the His41 catalytic (1.57% and 1.66%). Conversely, LMed-137 ( Figure 7B) performs several H-bonds, but all have a lifespan of less than 10%. Among them, three with the Thr26 residue (7.10%, 6.72%, and 4.29%), three with Asn51 (3.75%, 3.53%, and 2.59%), one with His41 (2.44%), three with Asn119 (2.83%, 2.47%, and 2.45%), two with Cys145 (1.44% and 1.27%), one with Tyr54 (1.09%), and three with Gln189 (1.55%, 1.04%, and 1.01%). Again, due to the short lifetime of these interactions, it was possible to confirm the detachment trend of these two ligands at the active site, as observed in the RMSD analysis.
Analyzing the ligand LMed-110 ( Figure 7C) showed two H-bonds interactions with a lifetime more significant (> 50%) with the residues His164 (69.39%) and Glu166 (57.48%). In addition, His164 performed another H-bonds interaction with the ligand's nitrogen (N1) atom with a short lifetime of 9.47%, and we also observed two H-bonds with a short lifetime of 4.51% and 1.23% with the residues His41 catalytic and Asn142, respectively.
The LMed-117 ( Figure 7D) performed H-bonds with only two residues. With the residue Asn142, the ligand performed one interaction of H-bond with a short lifetime (12.02%). After 27 ns, it is possible to observe one new H-bond with the catalytic histidine (His41) with a lifetime of 32.40%. Due  to the short lifetime of these interactions, it was possible to confirm the detachment trend of this ligand at the active site observed in the RMSD analysis. The protein-ligand binding free energy (DGbinding) can be defined as the summation of the energies associated with vacuum inter-action (vacuum potential energy) and solventdependent interactions (polar and nonpolar solvation energies). Therefore, DGbinding for each M pro -ligand complex was calculated considering the MMGBSA method (Kumari et al., 2014) for the most populated clusters (Table 2).
According to the DGbinding values, LMed-110 and LMed-136 would have better affinities for M pro than LMed-117 and the tautomer LMed-137. These results are similar to the docking study discussed above. Furthermore, for all inhibitors, the non-polar interactions (DEvdw þ DEsasa) would be more critical for the complex stability than electrostatic interactions (DEelect). The nonpolar and polar equilibrium could be related to the hydrogen bonding interactions since the high persistence of hydrogen bond interactions contribute more to ligand stability in the binding site (Albuquerque et al., 2020). The LMed-117 has a profile like the LMed-137; the low value of DGbbinging may be related to the lack of hydrogen bonding interactions with high persistence. Unfortunately, this behavior is not enough to keep the ligand stable at the active site. In general, this tendency to displace LMed-117 and LMed-137 structures from the active site is evidenced by the high value of RMSD. Therefore, it can be stated that due to the better binding affinity of the LMed-110 and LMed-136, both can act as a possible inhibitor for the M pro of SARS-CoV-2.

In silico drug-likeness, pharmacokinetic parameters (ADME) predictions
After the virtual screening based on molecular docking and MD simulations, we performed an in silico study of the physicochemical, pharmacokinetic, and toxicity properties to predict the behavior of these structures in living organisms (Table 3). Lipinski's rules (Lipinski et al., , 2001 are widely used to predict drug-likeness, i.e., the similarity of small molecules to known drugs related to absorption and oral bioavailability. After analyzing the physicochemical properties, LMed-136, LMed-137, LMed-110, and LMed-117 do not violate these rules, demonstrating that the compounds are similar to known drugs and could present good oral bioavailability such as Nirmatrelvir, a drug for oral use, which also does not violate Lipinski's rule. Besides, they fit Veber's parameters, indicating good membrane permeability (Veber et al., 2002). Furthermore, the LogS values for LMed-136, LMed-137, and LMed-110 were À 4.99, À 4.99, and À 4.96, respectively. They are similar to the drug Remdesivir (-4.12) with administration intravenously, showing that the compounds have a moderate solubility in water.
Regarding pharmacokinetic properties, our study revealed that the compounds have high gastrointestinal absorption, which indicates an ability to cross cell barriers in the GI tract, as does Nirmatrelvir. Concerning the blood-brain barrier, it allows the blood vessels that vascularize the central nervous system (CNS) to regulate the movement of ions rigidly, molecules, and cells between the blood and the brain, allowing an adequate neuronal function and protection of the neural tissue (Daneman & Prat, 2015). Therefore, analyzing our prediction, only LMed-110, and LMed-117 compounds can permeate the BBB, which can cause adverse or toxic effects in the CNS. P-glycoprotein (P-gp), an efflux transporter of toxins and xenobiotics out of cells, plays an essential role in the pharmacokinetics of several drugs, mainly interfering with drug absorption resulting in low bioavailability and/or nonlinear absorption (Ohashi et al., 2019). Therefore, no relevant changes in the absorption of these compounds are expected since they do not act as a substrate for P-gp. On the other hand, Cytochrome P450 (CYP) enzymes are essential in the metabolism and elimination of several xenobiotics, including drugs (Wasukan et al., 2019). Herein, we use a forecast model with five CYP450 isoforms of CYP1, CYP2, and CYP3 families (CYP1A2, CYP2C19, CYP2C9, CYP2D6, and CYP3A4). After these predictions, we found that LMed-136 and its tautomer LMed-137 are possible inhibitors of CYPC19 and CYP2C9 isoforms. On the other hand, LMed-110 and LMed-117 may have more significant drug interactions, as it probably inhibits four isoforms (CYP1A2, CYPC19, CYP2D6 and CYP3A4) similar to Remdesivir, which can inhibit three isoforms (CYP2C9, CYP2D6, and CYP3A4). Consequently, the inhibition of these isoforms can decrease the efficacy and increase the side effects and toxicity of co-administered drugs, indicating the need for further studies to verify an interaction with drugs metabolized by these CYP enzymes. In addition, we also verified whether the compounds LMed-136, LMed-137, LMed-110, and LMed-117 could act as a possible substrate for CYPs. This prediction is useful because if a substance behaves as a substrate for CYPs, it can lead to low oral bioavailability caused by pre-systemic metabolism. For example, Nirmatrelvir is a known substrate of CYP3A4, the most important drug-metabolizing enzyme (Klein & Zanger, 2013), which was corroborated by in silico predictions, so it must be administered in combination with Ritonavir to minimize its metabolism and increase its half-life (FDA, 2021c).
In contrast, our predictions showed that none of the LMed compounds are substrates of CYP3A4, which can be considered an advantage compared to Nirmatrelvir. However, our predictions showed that the compounds LMed-136, LMed-137, and LMed-110 could be possible substrates of the CYP2C9 and CYP2D6 isoforms, whereas LMed-177 is a possible substrate only of the CYP2D6 isoform. In addition, we predicted the positions of metabolically labile atoms in the substrate where metabolic reactions can be initiated by analyzing the sites of metabolism (SOMs). The study showed that the carbon of the guanidinium group is the most susceptible to oxidation by CYP2C9 and CYP2D6.

Conclusion
In the conclusion, our molecular docking-based virtual screening approach identified four benzoylguanidines (LMed-136, LMed-137, LMed-110, and LMed-117) with the highest binding affinity for M pro . This is the first report of benzoylguanidines as possible inhibitors against SARS-CoV-2 M pro . Furthermore, molecular dynamics simulations (150 ns) of the four main anchored structures revealed essential insights into the possible molecular interactions with the molecular target, such as the conformational change concerning molecular docking results and the short duration of H-bonds. After these analyses, LMed-110 and LMed-136 were elected as potential hits, considering that they performed H-bonds with a lifetime more significant than 50%, which indicates a greater affinity for the active site. Therefore, we evaluated the pharmacokinetic properties of our compounds using ADME in silico predictors, which indicate that LMed-136 can be used for medicinal purposes due to its drug-like properties, good pharmacokinetics. Our results point to LMed-136 as a promising compound, guiding us towards the following steps to be carried out in this work, such as in vitro cytotoxicity and antiviral assays, as well as the design of new possible inhibitors by medicinal chemistry strategies.

Acknowledgments
The Molecular Dynamics simulations were realized using resources of the "Centro Nacional de Processamento de Alto Desempenho em São Paulo (CENAPAD-SP)."

Disclosure statement
No potential conflict of interest was reported by the author(s).