In Silico Study of Coumarins and Quinolines Derivatives as Potent Inhibitors of SARS-CoV-2 Main Protease

The pandemic that started in Wuhan (China) in 2019 has caused a large number of deaths, and infected people around the world due to the absence of effective therapy against coronavirus 2 of the severe acute respiratory syndrome (SARS-CoV-2). Viral maturation requires the activity of the main viral protease (Mpro), so its inhibition stops the progress of the disease. To evaluate possible inhibitors, a computational model of the SARS-CoV-2 enzyme Mpro was constructed in complex with 26 synthetic ligands derived from coumarins and quinolines. Analysis of simulations of molecular dynamics and molecular docking of the models show a high affinity for the enzyme (∆E binding between −5.1 and 7.1 kcal mol−1). The six compounds with the highest affinity show K d between 6.26 × 10–6 and 17.2 × 10–6, with binding affinity between −20 and −25 kcal mol−1, with ligand efficiency less than 0.3 associated with possible inhibitory candidates. In addition to the high affinity of these compounds for SARS-CoV-2 Mpro, low toxicity is expected considering the Lipinski, Veber and Pfizer rules. Therefore, this novel study provides candidate inhibitors that would allow experimental studies which can lead to the development of new treatments for SARS-CoV-2.


INTRODUCTION
In recent years, different viruses have emerged in around of the world. These diseases are a generation of respiratory diseases in infected patients, also due to the rapid dissemination of the diseases. These kinds of viruses include the severe acute respiratory syndrome coronavirus (SARS-CoV), Middle East respiratory syndrome coronavirus (MERS-CoV), avian influenza A/H 7 N 9 and H 5 N 1 viruses, and Nipah virus (Yuen et al., 1998;Peiris et al., 2003;MacNeil and Rollin, 2012;Marsh and Wang, 2012;To et al., 2012;To et al., 2013;Zaki et al., 2012;Chan et al., 2015). The capacity of these viruses to evolve and infect humans has been associated with the close interaction occurring between human populations and different animal species in markets of densely populated areas (Chan et al., 2015). In December 2019, cases of atypical pneumonia began to be observed in the city of Wuhan (China) (Lu et al., 2020;Zhu et al., 2020). By January 2020, the etiological agent was classified as a new member the of family Coronaviridae and genus β-coronavirus (2019-nCoV) that differ from SARS-CoV and MERS-CoV. The genome of 2019-nCoV shares an 82% sequence identify to SARS-CoV (Elfiky, 2020;Hui et al., 2020;Rothan and Byrareddy, 2020;World Health Organization, 2020). As a matter of fact, its genome has high similarity with the genome of a bat coronavirus (96.2% identity), which has allowed the virus to be associated with a zoonotic origin (Lu et al., 2020;Wu et al., 2020;Zhou et al., 2020). According to the International Committee Virus Taxonomy the new β-coronavirus was called severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and on February 11 the set of symptoms associated with this new virus was designated as COVID-19 by the World Health Organization (WHO) . The evolution of infections at a global level increased as the number of cases and deaths, with the most affected countries being the USA, India, Brazil, Russia, France and United Kingdom (COVID-19 Map -Johns Hopkins Coronavirus Resource Center, 2020), which together have presented more than 50% of the global cases in more than 180 countries, declaring it a pandemic on March 11 (WHO Director-General's opening remarks at the media briefing on COVID-19-20 March, 2020).
Like all other coronaviruses, SARS-CoV-2 is composed by single-stranded RNA as their genetic material with an approximate length of 29,891 nucleotides and a 5′-cap structure and 3′-poly-A tail, encoding 9,860 amino acids (Chan et al., 2020). This RNA encodes both the structural and non-structural proteins of the virus. Among the structural proteins, there is the Spike (S) (present in all coronaviruses), Nucleocapsid (N), Matrix (M), and Envelope (E) (Chan et al., 2020;Walls et al., 2020;Wrapp et al., 2020). Proteases and RNAdependent RNA polymerase constitute the non-structural proteins of the virus. The genome contains at least six open reading frames (ORFs), the first of these ORF occupies about 60% of the length of the genome and translates two polyproteins known as pp1a and pp1ab, which are processed by the main protease (M pro , also called 3CL pro ) and papain-like proteases (PLPs) (Jin et al., 2020a;Zhang et al., 2020b;Chen et al., 2020). Consequently, inhibiting M pro activity blocks virus replication and thus affects the life cycle of SARS-CoV-2. Compounds derived from coumarins and quinoline have been tested against various viruses (McKee et al., 2020;Mishra et al., 2020). Quinolines have recently been used in experimental treatments for SARS-CoV-2 infected persons in several countries (Arshad et al., 2020;Lopez et al., 2020;Mori et al., 2020) and it has been proposed that coumarins inhibit the replication of several viruses including influenza (Pavurala et al., 2018), HIV (Jesumoroti et al., 2019), Dengue (Coulerie et al., 2013), Chikungunya (Hwu et al., 2019), hepatitis (Hwu et al., 2008) and filoviruses (EBOLA, Marburgvirus (MARV) and Cuevavirus) (Liu et al., 2019).
In this study, we evaluated twenty-sixth molecules derived from coumarins and quinolines as promising SARS-CoV-2 M pro inhibitors, and so, by using computational biochemistry protocols we tried to find the most appropriate molecules that can act as potential anti-SARS-CoV-2 activity drugs. Six of the compounds evaluated are highlighted, which are CTR9, 7HC6, CTR6, 7HC5, 7HC3 and 8HQ6. We performed molecular docking (rigid), efficiency calculations of ligands, pharmacological and toxicological property predictions (ADMET), and molecular dynamics simulations (MD) simulations, together with MM-GBSA binding free energy predictions to identify the binding characteristics for identifying the inhibitors of SARS-CoV-2 M pro .

Compounds Set
In this work, we use twenty-sixth ligands selected for their possible capability to inhibit SARS-CoV-2 protease M pro . The coumarins and quinolines derivatives were extracted from two compound series that have been synthesized by the laboratories of applied chemistry of the Universidad de Ibagué, in Ibague-Colombia. The work was based on strategies that the authors have reported in the literature (García-Beltrán et al., 2013;Mena et al., 2015;Aguirre et al., 2017;Garćia-Beltran et al., 2017). 13b is our reference molecule and was obtained from the Protein Data Bank (PDB) ) (PDB id: 6Y2F). These molecules were designed in silico and evaluated using docking methodologies and physicochemical and pharmacokinetic descriptors, and to predict ADME parameters. Their chemical structures are shown in Figure 1. Their molecular conformations were optimized using PM6-D3H4 semi-empirical method (Stewart, 2007;Řezáč and Hobza, 2012) as implemented in MOPAC2016 (Stewart, 2016) software. The optimized molecules were used for molecular docking simulations in order to study the interactions established by these compounds in the SARS-CoV-2 M pro pocket.

Molecular Docking
We made molecular docking analyses to examine the potential binding modes of ligands to the main protease M pro of SARS-CoV-2, as potential inhibitors. Then, based on structural information obtained from the crystal structures of M pro in complex with other ligands we established the binding site of the proposed inhibitors for SARS-CoV-2 M pro . In addition to delimit the binding region of possible inhibitors to M pro (Dai et al., 2020;Jin et al., 2020a;ul Qamar, 2020;Xu et al., 2020;Zhang et al., 2020a. We established the key residues for catalysis according to experimental and theoretical data. Finally, we used AutoDock (v 4.2.1) and AutoDock Vina (Trott and Olson, 2010) for all dockings in this study. The initial 3D inhibitors structures were drawn using Discovery Studio (Dassault Systèmes BIOVIA, 2017) 3.1 (Accelrys, CA) which were optimized (considering the RMS gradient of 0.001 kcal/mol) using the PM6-D3H4 semi-empirical method (Stewart, 2007;Řezáč and Hobza, 2012) implemented in the MOPAC2016 (Stewart, 2016) software. PM6-D3H4 introduces dispersion and hydrogen-bonded corrections to the PM6 method. The ligand files were prepared using the AutoDockTools package (Sanner, 1999) provided by AutoDock through accepting all rotatable bonds; moreover, the atomic charges are computed toward the PM6-D3H4 procedure, and non-polar hydrogen atoms are merged. The semi-empirical method has shown to increase significantly Frontiers in Chemistry | www.frontiersin.org February 2021 | Volume 8 | Article 595097 3 docking accuracy and cluster population of the most accurate docking (Bikadi and Hazai, 2009;Fanfrlík et al., 2010;Hou et al., 2013). The crystal structure of SARS-CoV-2 M pro (PDB Code: 6YB7, https://www.rcsb.org/structure/6YB7), was downloaded from the PDB (Berman et al., 2000). The resolution of the retrieved structure was 1.25 Å. The SARS-CoV-2 M pro was treated with the Schrödinger's (Schrödinger, 2020) Protein Preparation Wizard (Madhavi Sastry et al., 2013); polar hydrogen atoms were added, non-polar hydrogen atoms were merged, and charges were assigned. Docking was treated as rigid and carried out using the empirical free energy function and the Lamarckian Genetic Algorithm provided by AutoDock Vina (Morris et al., 1998). The docking grid dimensions were 30.75 × 30.75 × 30.75 Å, making the binding pocket of SARS-CoV-2 M pro the center of mass between amino-acid residues (Cys145 and His41) of the catalytic site. All other parameters were set as defined by default through AutoDock Vina. Dockings were repeated 50 times with space search exhaustiveness set to 100. The best interaction binding energy (kcal·mol −1 ) was selected for evaluation and analyzed according to the potential intermolecular interactions (protein/ligand), such as hydrogen bonding, hydrophobic interactions and the cation-π, π-π stacking. Docking results 3D representations were used. VMD molecular graphics system (Humphrey et al., 1996).

Ligand Efficiency Approach
Ligand efficiency (LE) calculations were performed using one parameter K d . The K d parameter corresponds to the dissociation constant between a ligand/protein, and their value indicates the bond strength between the ligand/ protein (Abad-Zapatero, 2007, Abad-Zapatero, 2013Abad-Zapatero et al., 2010). Low values indicate strong binding of the molecule to the protein. K d calculations were done using the following equations: where ΔG 0 corresponds to binding energy (kcal mol −1 ) obtained from docking experiments, R is the gas constant whose value is 1.987207 cal/mol K and T is the temperature in degrees Kelvin. At standard conditions of aqueous solution at 298.15 K, neutral pH and remaining concentrations of 1 M. The ligand efficiency (LE) allows us to compare molecules according to their average binding energy (Reynolds et al., 2008;Abad-Zapatero, 2013). Thus, it determined as the ratio of binding energy per nonhydrogen atom, as follows (Abad-Zapatero, 2007, Abad-Zapatero, 2013Abad-Zapatero et al., 2010;Cavalluzzi et al., 2017): where K d is obtained from eq. (2) and HAC denotes the heavy atom count (i.e., number of non-hydrogen atoms) in a ligand.

ADMET Properties
The purpose of calculating ADMET profiles is to supply, with reasonable accuracy, a preliminary prediction of the in vivo behavior of a compound to assess its potential to become a drug (Yu and Adedoyin, 2003). The molecules used in this study were submitted to the calculation of their absorption, distribution, metabolism, excretion and toxicological properties (ADMET). Also, the physicochemical properties such as molecular hydrogen bond acceptor (HBA), hydrogen bond donor (HBD), molecular weight (MW), topological polar surface area (TPSA), rotatable bond count (RB) and octanol/ water partition coefficient (LogP) were calculated using SwissADME webserver (Daina et al., 2017). Compound toxicological properties were analyzed taking into account the Lipinski, Veber and Pfizer toxicity empirical rules, see Table 1.

Molecular Dynamics Simulations
MDs calculations were performed for the lowest six binding energy docking and the compound 13b, which is our reference ligand. These calculations were also obtained from Protein Data Bank (PDB id: 6Y2E). The ligands were bound to SARS-CoV-2 M pro protein (PDB ID:6YB7) in aqueous solutions with an explicit solvent TIP3P water model (Neria et al., 1996) (≈16.000 water molecules). Protonation states of ionizable residues corresponding to pH 7.0 were determined by H++ web interface for computes pK values of ionizable groups in macromolecules and adds missing hydrogen atoms according to the specified pH of the environment (Anandakrishnan et al., 2012). Besides, NaCl ions were modeled to neutralize the systems and maintain an ionic concentration of 0.15 mol/L. The compounds were parameterized by GAFF Force Field for organic molecules. (Wang et al., 2004;ÖzpInar et al., 2010), using the Antechamber module in AmberTools18 with AM1-BCC charges, (Jakalian et al., 2000). The protein structures were modeled with the force field ff14SB. (Salomon-Ferrer et al., 2013). The simulations were carried out using a standard MD protocol: (I) Minimization and structural relaxation of water molecules with 2000 steps of minimization (downward step) and MD simulation with an NPT (300 K) assembly by 1,000 ps using harmonic restrictions of 10 kcal molÅ −2 for protein and ligand; (II) minimization of the complete  structure considering 2000 downstream minimization steps and 6,500 steps of conjugate gradient minimization; (III) the minimized systems were progressively heated to 300 K, with harmonic restrictions of 10 kcal mol Å −2 in the carbon skeleton and ligand during 0.5 ns; (IV) the system was then balanced by 0.5 ns maintaining the restrictions and then by 5 ns without restrictions to 300 K in a canonical assembly (NVT); and (5) finally, a production dynamic was carried out with an isothermal isobaric assembly (NPT) without restrictions for 200 ns at 310 K and 1 atm with a temporary passage of 2 fs. In the MD simulation, the temperature was controlled by the Langevin dynamics with a collision frequency of 1 ps −1 (NVT) and the pressure with the Berendsen barostat (NPT). Besides, the Particle Mesh Ewald (PME) method with a cut-off value of 10 Å was used to treat nonbonding and long-range electrostatic interactions. All MD simulation calculations were performed using the AMBER-GPU Implementations18 (Mermelstein et al., 2018). Molecular visualization of the systems and MD trajectory analysis was carried out with the VMD software package (Humphrey et al., 1996).

Free Energy Calculation
The molecular MM/GBSA method was employed to estimate the binding free energy of the protease-ligand complexes. For calculations from a total of 200 ns of MD, the last 50 ns were extracted for analysis, and the explicit water molecules and ions were removed. The MM/GBSA analysis was performed on three subsets of each system: the protein alone, the ligand alone, and the complex (protein-ligand). For each of these subsets, the total free energy (ΔG tot ) was calculated as follows: where H MM is the bonded and Lennard-Jones energy terms; G solv is the polar contribution of solvation energy and nonpolar contribution to the solvation energy; T is the temperature; and ΔS conf corresponds to the conformational entropy (Hayes and Archontis, 2012). Both H MM and G solv were calculated using AMBER 18 program with the generalized Born implicit solvent model (Götz et al., 2012;Song et al., 2019). ΔG tot was calculated as a linear function of the solvent-accessible surface area, which was calculated with a probe radius of 1.4 Å (Abroshan et al., 2010). The binding free energy of SARS-CoV-2 M pro and ligand complexes (ΔG bind ) were calculated by the difference where G_tot values are the averages over the simulation.

Non-Covalent Interactions
The principal cluster of main component analysis of trajectory were analyzed with the non-covalent interaction index (NCI) (Johnson et al., 2010;Contreras-García et al., 2011) using NCIPLOT program (Contreras-García et al., 2011 to identify and map non-covalent interactions, such hydrogen bonds, steric repulsion, and van der Waals interactions, using the promolecular densities (ρ pro ), computed as the sum of all atomic contributions. The NCI is based on the electron density (ρ), its derivatives and the reduced density gradient (s). The reduced density gradient is given by: These interactions are local and manifest in real space as low-gradient isosurfaces with low densities which are interpreted and colored according to the corresponding values of sign(λ 2 )ρ. The surfaces are colored on a bluegreen-red scale according to the strength and type of interaction. Blue indicates strong attractive interactions, green indicates weak van der Waals interactions, and red indicates a strong non bonded overlap.

Molecular Docking Analysis
The docking results, which were conducted to estimate the possible binding sites of potential inhibitors on SARS-CoV-2 M pro . The genetic material of SARS-CoV-2 expresses multiple proteins (more than 20 proteins), among these proteins the main protease (Mpro) is identified, a molecule similar to 3 chymotrypsin (3CL pro ) that shows a similarity of 96.1% with the 3CL pro of SARS-CoV. 3CL pro plays a very important role in replication and transcription processes of the virus genome (Hui et al., 2020). Therefore, 3CL pro is a strategic drug target in the inhibition of the SARS-CoV cycle. The protease is active as a homodimer, structured by the dimerization of two protomers designated as monomer A and monomer B, and the catalytic dyad in each protomer is defined by Cys145 and its residues . This has led to the development of multiple studies with experimental and computational approaches in search of possible inhibitors that can effectively block the activity of this protease (Balaramnavar et al., 2020;Dai et al., 2020;Gao et al., 2020;Jin et al., 2020a;Jin et al., 2020b;Ngo et al., 2020;Zhang et al., 2020b). Work with the 3C-like proteinase from SARS coronavirus revealed that the Cys145 residue is key at the active site of 3CLpro (Huang et al.,  2004), this advance allowed the mentioned residue to be an attractive target for covalent ligands to bind acting as inhibitors of 3CL pro . This amino acidic residue is also a popular target for covalent inhibitors because of its intrinsic reactivity at physiological pH (Cuesta et al., 2020). Tung Ngo et al. report in recent studies that additionally Glu166 residue has a prominent and important role in binding ligands to SARS-CoV-2 M pro (Ngo et al., 2020). Twenty-seven inhibitors, including the reference ligand 13b were evaluated in silico anti-SARS-CoV-2 activity. The results of this study of molecular docking calculations indicate the strong interactions of molecules derived from coumarins and quinolones targets the Cys-His catalytic dyad (Cys145 and His41) in the binding pocket of SARS-CoV-2 M pro . The results of the binding of these molecules are presented in Table 2. Meanwhile, the binding to various amino acid residues due to their presence in the conserved region of the active site in all compounds is seen and presents a very important role in enzymatic catalysis.
Docking results with SARS-CoV-2 M pro , indicated that all ligands present binding energies between −5.3 and −7.2 kcal mol −1 ( Table 2), with a difference of 0.1 kcal mol −1 between 13b and CTR6, which means that the level of stability is very similar between these two protease-complex ligands. AutoDock Vina, presents that the ligand 13b attached to the co-crystal, and the re-coupled ligand 13b presents an RMSD value of 3.1 Å, suggesting a partially acceptable value of the coupling method. Furthermore, it shows that the compound CTR6 gives the lowest energy (−7.1 kcal-mol −1 ) in complex with the protease, which is the best score when compared to other docked compounds used in this study. CTR6 gives better score than 7HC6 (−6.7 kcal mol −1 ), 7HC5 (−6.6 kcal mol −1 ), CTR9 (−6.6 kcal mol −1 ), 7HC3 (−6.6 kcal mol −1 ), 8HQ6 (−6.5 kcal mol −1 ) and the other compounds. The coumarins and quinolines derivatives are located inside the protein pocket, by means of electrostatic and hydrophobic interactions with the residues Arg188, Asn142, Asp187, Cys145, Gln189, Glu166, Gly143, His164, His163, His41, Leu167, Leu27, Met165, Met49, Phe140, Pro168, Ser144, Thr25 and Thr26. In relation to 13b, it shows two hydrogen bonding interactions with Leu167 and Glu166 residues. In the case of CTR6, 7HC6 and CTR9 in complex with SARS-CoV-2 M pro (Figure 2 and Table 2) show possible hydrogen-bonding interactions with the residues Ser46, Met49, Asn142, Met165, Glu166, . This allows us to conclude that the hydrogen bonds and hydrophobic forces, are the majority interactions that dominate these complexes. Interactions between the rest of the compounds and SARS-CoV-2 M pro are reported in Table 2.

Ligand Efficiency Analysis
The parameters dissociation constant (K d ) and ligand efficiency (LE) were used to compare the affinity of the molecules studied

Molecular Dynamics Simulation and MM/ GBSA Analysis
Molecular dynamics simulations were performed in 200 ns to analyze the steady nature and conformations stability of ligand-SARS-CoV-2 M pro complexes (ligands: CTR9, 7HC6, CTR6, 7HC3, 8HQ6 and 7HC5). The RMSD was used to estimate the stability of protein-ligand systems. RMSD trajectories of SARS-CoV-2 M pro -ligand complexes during 200 ns simulation indicated that the complexes formed with the ligands during the DM simulations have a high stability during the simulation time (Figure 3). The structure does not show significant changes, in this case there is an increase of the RMSD until reaching a point in which the values fluctuate around values of 0.5 and 1.8 Å of RMSD. After a Molecular dynamics of 200 ns the structures remain within the parameter that considers the system to be in equilibrium, therefore, no complex suffered structural destabilization during the simulation. Deviations with a maximum difference of 3.0 Å of RMSD (Carugo, 2003) indicate that the system is in equilibrium, situation that is fulfilled for the simulation of the possible SARS-CoV-2 M pro indicating equilibrium states of the ligands within the active site of protease. Also, RMSD curves for 7HC5 and CTR6 are remarkably more stable than those of CTR9, 7HC6, 7HC3 and 8HQ6. To complement the analysis carried out from the RMSD, the study of the Radius of Gyration (R Gyr ) was carried out for the same runs. From R Gyr analysis Figure 4, we can conclude that the R Gyr of ligands 7HC5 and 7HC3 have values that oscillate in an interval close to 3.0 Å, for the case of CTR9, 8HQ6 and 7HC6 have values fluctuating in an interval close to 3.5 Å and for the CTR6 ligand it has values higher than 4.0 Å. The stable values during the 200 ns simulation for R Gyr indicate that ligand binding at the active site of the protein does not induce major conformational changes in the protein structure.
With the purpose of identifying the deviation of the ligand from respect to its initial position and the movement of proteins residues, the Normal Mode Analysis (NMA) and Root Mean Square Fluctuations (RMSF) values were calculated averaging over all the conformations sampled during the last 100 ns simulation. The NMA and RMSF were calculated using the Cα atom of each amino acid residue as a reference and the graph was used to represent the fluctuations in the residue level. NMA plot in Figure 5 shows a similar trend of residue fluctuation profile for the complex with an average NMA of 1.0 × 10 -6 fluctuations. The two CTR6-protease and 7HC3-protease complexes showed a comparatively higher fluctuation in some residues. This trend in the quadratic displacement figure of the complex suggests that the binding of the six compounds to the protein showed stability and no effect on the flexibility of the protein was observed in the whole range of the simulations. The N-terminus to C-terminus of the proteins normally present great fluctuations and in most cases their movement does not represent importance. As shown in Figure 5, the RMSF graph, all the compounds made the residues present the same fluctuation, except for the two 8HQ6-protease and 7HC5-protease complexes showed a comparatively higher fluctuation in some residues. We consider and from the theoretical point of view, that these ligands present movements in the active site to settle in the best orientation, reason why part of the amino acidic residues fluctuate more than normal. However, other ligands exhibit good behavior and their RMSF values show that they can handle the fluctuations of the residues. On the active site of the protease, the fluctuation values of the main residues (Zhang et al., 2020a) (His41, His163, His164, Phe140 and Cys145) of the six selected molecules were similar among them. The results of the MD simulation indicate that two of the ligands obtained from the coupling analysis (CTR9 and 7HC6) remain close to their initial locations even in uncontrolled simulations, which points to the constitution of stable complexes. From these results it can be clearly deduced that it is likely that the molecules CTR9 and 7HC6 play the same role in inhibiting SARS-CoV-2 M pro as 13b.
The analyses of trajectories indicate that during most of the simulation the ligands CTR9, 7HC6, CTR6, 7HC5, 7HC3 and 8HQ6 maintain hydrogen bonds with residues of the active site of SARS-CoV-2 M pro . However, the number of hydrogen bridges formed was different for each ligand (Figure 6). CTR9 formed three hydrogen bridges between the residues Cys44, His41 and Thr26, highlighting the participation of the residues Glu166, Cys145, Asn142 and Asp187. 7HC6 formed two hydrogen bridges between the residues Thr24 and Thr25, highlighting the participation of the residues Leu27, Ser46, Met49 and Glu166. In the case of CTR6, three hydrogen bridges between the residues Glu166, Gln189 and Thr190 were determined, highlighting the participation of the residues Arg188 and Ser46. 7HC3 formed one hydrogen bridge with the residue Glu166, highlighting the participation of the residues His163, Met49 and Phe140. Three hydrogen bridges are formed between 8HQ6 and the residues Asp187, Glu166 and His163, highlighting the participation of the residues Gln189, Phe140, Ser144 and Val186. Finally, 7HC5 formed one hydrogen bridge with the Glu166 residue, highlighting the participation of the residues Leu50, Asn142 and Phe140. These residues, see Figures 6 and 7, are consistent with previous theoretical-experimental studies carried out by Dai et al. (Dai et al., 2020), where they detail the interaction that some of the synthesized compounds have with the active site of the protease.
The noncovalent interactions analysis labeled all the hydrogen-bonding interactions in total agreement with the molecular dynamics simulations, providing a qualitative confirmation of these interactions, using a topological and  visual analysis of a scalar field related to the electron density ( Figure 8). These results suggest that the better inhibitors character is due to direct mechanisms. Finally, the binding free energy (MM-GBSA) was estimated subsequent to the MD simulation; the last 50 ns for all the complexes and the results are given in Table 3. CTR9 and 7HC6 compounds depicted the lowest binding free energy (−24.8 and −24.7 kcal mol −1 ) with SARS-CoV-2 M pro , while the compounds CTR6, 7HC5, 7HC3 and 8HQ6 showed relatively higher binding energy (−24.7, −22.8, −22.6, −20.5 and −20.5 kcal mol −1 ). On the other hand, the reference compound 13b showed the lowest binding free energy (−29.1 kcal mol −1 ) with the SARS-CoV-2 M pro in comparison to compound CTR9, with a slight difference of 4.3 kcal mol −1 . Although compound 13b has a lowest binding free energy, it presents the problem of breaking the rules of Lipinski and Veber rules, however compound CTR9 does not.

ADMET Properties
In the search for new drugs, safety is very important and the regulations related to ADMET (Absorption, Distribution, Metabolism, Excretion and Toxicity), most of the time are the cause for a drug to fail. Therefore, it is of utmost importance to identify aspects such as the toxicity of compounds in early stages of development and thus avoid the loss of resources and time. (Sivamani et al., 2012). To evaluate the best ligands as potential anti-SARS-CoV-2 activity drugs; we have calculated some pharmacokinetic properties (Table 4). These results were compared to Veber's (Veber et al., 2002), Pfizer's (Hughes et al., 2008) and Lipinski's rule (Lipinski et al., 2001). in the development of new drugs if a molecule complies only with one of the Lipinski's rule is not an appropriate candidate and it is not relevant to continue with the study, however, by presenting a greater number of rules the probabilities of being a candidate begin to increase and deepen their study. In accordance with Veber's rule, if a compound does not satisfy at least two parameters, it is not a candidate for further development. In addition, Pfizer 3/75 toxicity rules have also been taken into account in this study, concluding that if any of the proposed ligands do not meet the established parameters, then it is not a suitable candidate. ADME prediction showed that in most cases, all the compounds proposed in this study satisfy with the Veber's, Pfizer's and Lipinski's rule. This suggests that these ligands could be safe molecules for use as anti-SARS-CoV-2 activity drugs. In the case of reference ligand 13b presents a violation of Lipinski's and Veber's rule, due to their molecular weight, topological polar surface area and rotatable bond count. The values of these properties are higher than the admissible limit, making this substance fat-soluble which indicates a tendency to be more toxic and less selective to their target. In the case of the six best compounds found in the docking simulations (CTR6, 7HC6, 7HC5, CTR9, 7HC3 and 8HQ6), all of them satisfactorily meet the Veber's, Pfizer's and Lipinski's rule. These compounds represent most promising compounds to molecular dynamics simulation and MM-GBSA.

CONCLUSION
This paper predicts that compounds derived from coumarins and quinolines that can be successfully potential drugs to treat viral diseases such as COVID-19. Herein, we used a computational chemistry protocol to identify the ligands most promising candidates that may inhibit main protease of SARS-CoV-2 activity determined by means of this protocol involves of molecular dockings, molecular dynamics simulations, MM-GBSA, NCI and ADMET properties to predict whether these compounds are appropriate to be utilized in an anti-COVID-19 therapy. We identified six compounds (CTR9, 7HC6, CTR6, 7HC5, 7HC3, 8HQ6) that are already synthesized (García-Beltrán et al., 2013;Mena et al., 2015;Aguirre et al., 2017;Garćia-Beltran et al., 2017) with a potential inhibition of main protease of SARS-CoV-2. These compounds might be repurposed against COVID-19. These hits were described as drug-like compounds and showed harmless ADMET properties and may aid in developing and optimizing more efficient and potent COVID-19 inhibitors. Trajectory analysis showed that the studied complexes display structural stability during the MD runs. These results encourage further in vitro and in vivo investigations and also preventively boost the traditional use of coumarins and quinolines derivatives preventively. We anticipate that the insights obtained in the present study may prove valuable for researching and developing novel anti-COVID-19 therapeutic agents in the future.