Discrimination of Naturally-Occurring 2-Arylbenzofurans as Cyclooxygenase-2 Inhibitors: Insights into the Binding Mode and Enzymatic Inhibitory Activity

2-arylbenzofuran-containing compounds are chemical entities that can be naturally produced by several organisms. A wide-range of activities is described for several compounds of this kind and they are, therefore, valuable moieties for a lead finding from nature. Although there are in-vitro data about the activity of 2-arylbenzofuran-related compounds against cyclooxygenase (COX) enzymes, the molecular level of these COX-inhibiting constituents had not been deeply explored. Thus, 58 2-arylbenzofurans were initially screened through molecular docking within the active site of nine COX-2 crystal structures. The resulting docking scores were statistically analyzed and good reproducibility and convergence were found to discriminate the best-docked compounds. Discriminated compounds exhibited the best performance in molecular dynamics simulations as well as the most-favorable binding energies and the lowest in-vitro IC50 values for COX-2 inhibition. A three-dimensional quantitative activity-structure relationship (3D-QSAR) was also demonstrated, which showed some crucial structural requirements for enhanced enzyme inhibition. Therefore, four hits are proposed as lead structures for the development of COX-2 inhibitors based on 2-arylbenzofurans in further studies.


Introduction
Cyclooxygenase, COX, (or Prostaglandin H2 synthase), is a well-known enzyme that catalyzes the rate-limiting step in the formation of prostanoids such as prostaglandins (PGs) and thromboxane A2 (TxA2). These products are considered as important mediators of fever, pain, and inflammation [1]. Three COX isozymes are known such as the constitutive form including COX-1 (which exists in several tissues), the inducible form, COX-2 (which is expressed throughout inflammation), and the COX-1 splice variant, COX-3 [1]. In recent decades, the fact for understanding the distinct functions of COX isozymes has encouraged many studies, including those dedicated to the search for anti-inflammatory agents. In this sense, it is well-recognized that selective COX-2 inhibitors have been developed as securer substitutes for the nonsteroidal anti-inflammatory drugs (NSAIDs), mainly due to the adverse effects caused by the inhibition of the COX-1 enzyme. Therefore, a pronounced interest for finding selective anti-COX-2 agents still remains [2].
Several types of compounds (natural or synthetic) with anti-COX-2 properties have been discovered in recent years. Several chemical entities have been studied as selective COX-2 inhibitors with interesting and important potential [3]. However, the most important ones are those named coxibs (such as celecoxib, rofecoxib, and valdecoxib), which have been approved due to the lower adverse effects in comparison to the traditional NSAIDS (non-steroidal anti-inflammatory drugs) [3]. However, it has recently been reported that some coxibs and other selective COX-2 inhibitors cause non-favorable gastrointestinal and cardiovascular events [3,4]. Because of that, novel agents for COX-2 inhibition are still required.
A huge interest on benzofuran-containing compounds has been originated several years ago because of their important biological and pharmacological properties and their extensive natural occurrence [5]. Benzofuran moiety is considered to be an indole isostere (a potent pharmacodynamic moiety exhibiting anti-inflammatory activity) and, consequently, several anti-inflammatory benzofuran-based compounds have been studied [6,7]. However, despite the plethora of reported benzofuran derivatives, 2-arylbenzofurans have not been sufficiently explored at a molecular level for COX-2 inhibition in spite of some of them exhibiting selective in-vitro anti-COX-2 activity [8,9]. Thus, as part of the research on anti-inflammatory compounds of natural origin, a combined in-silico study comprising different computational approaches/tools were used on the basis of three motivations: (1) discriminating putative COX-2 binders through a molecular docking-based screening of 2-arylbenzofuran-related compounds, within the active site of nine reported crystal structures of COX-2, (2) describing some insights into the binding mode of the best-docked compounds using molecular dynamics and binding energy calculations, and (3) recognizing some structural requirements of the best-docked compounds toward COX-2 inhibition by means of comparative molecular field analysis (CoMFA). This last purpose required the expansion of the in-silico predictions to experimental results through the in-vitro COX-2 inhibitory assay performed on an in-house collection of 2-arylbenzofurans. Results revealed four very important hits to be used in future studies for developing novel selective COX-2 inhibitors based on 2-arylbenzofurans.

Protein Preparation
The X-ray crystallographic structure of nine COX-2 proteins (from Mus musculus) were obtained from the Protein Data Bank (PDB) (http://www.rcsb.org) at a resolution < 3.00 Å (Table 1). Water molecules, ligands, and other heteroatoms were removed from the protein molecules. The addition of hydrogen atoms to the proteins was then performed. These crystallographic structures were retained without any processing for molecular docking. The co-crystallized inhibitors/substrates were employed to define the corresponding active site, the flexible residues within this active site, and a validation criterion of docking calculations (re-docking). By means of a literature survey, a maximum of five common residues involved in the binding site were defined as flexible residues: Arg106, Val335, Ser339, Tyr341, and Ser530. PDB files of proteins were effectively prepared using the AutoDock/Vina (Molecular Graphics Lab, The Scripps Research Institute, La Jolla, CA, USA) plugin under PyMOL 2.3 (Schrödinger LLC, New York, NY, USA), and saved as pdbqt files. COX-2 structures were also superimposed by sequence alignment for calculating the root-mean-square deviation (RMSD) and sequence similarity using Discovery Studio Client 16.1 (Biovia, San Diego, CA, USA).

Ligand Preparation
A targeted, custom-made library comprising only 2-arylbenzofurans (n = 58) was compiled from the information of an in-house collection of isolated compounds as well as different published works reporting anti-inflammatory activity of structurally-related 2-arylbenzofurans of natural origin. Compiled compounds were three-dimensionally (3D)-sketched using the molecular builder and visualization tool Avogadro 1.2.0 (open-source, http://avogadro.cc). A Monte-Carlo randomized conformational search was then performed, without any geometrical restrictions, using the semi-empirical AM1 parametrization implemented in the Spartan'14 software (Wavefunction, Irvine, CA, USA) with a limit of 500 conformers. Energetically-lowest stable conformers, within a 6 kcal/mol energy range, were geometrically optimized at a density functional theory (DFT) level using B3LYP hybrid functional and 6-311++G (2df, 2p) basis set that is also implemented in Spartan'14. After structural optimization, harmonic vibrational frequencies were calculated at the same level of theory to verify the reliability of the stationary point at a minimum and all structures converged successfully. During geometry optimization, the highest energy occupied molecular orbital (HOMO) and lowest energy unoccupied molecular orbital (LUMO) were also computed for further analysis. The verified, optimized ligand structures were saved as PDB-files for further molecular docking and molecular dynamics analysis. Additionally, structural similarity analysis of 2-arylbenzofurans was performed using the substructure fragment dictionary-based binary fingerprint descriptor (FragFp) implemented in the open-source data visualization and analysis program DataWarrior [17]. The absorption, distribution, metabolism, and excretion-toxicity (ADMET) properties of test compounds were computationally predicted using the Free ADMET Filtering tool (FAF-Drugs4, http://fafdrugs4.mti.univ-paris-diderot.fr).

Molecular Docking
Molecular docking calculations were performed by means of Autodock/Vina (Molecular Graphics Lab, The Scripps Research Institute, La Jolla, CA, USA) using PyMOL 2.3 (Schrödinger LLC, New York, NY, USA) and running on a Ubuntu 12.04 server equipped with a dual Intel Xeon ® processor CPU @ 2.6 GHz (32 CPU), 64 GB DDR3 RAM. Prepared proteins and ligands were used for simulations using flexible residues. The selection of such flexible residues within the active site was based on a position at 4.0 Å from co-crystallized ligand (i.e., substrate or inhibitor, as presented in Table 1 and Figure A1). Thus, the active site was placed into a 24 × 24 × 24 Å cubic box at the geometrical center of the set of selected flexible residues, which has 0.375 Å as grid point spacing. Ten different pdbqt-files of each ligand were used as replicates to be docked into the flexible active site of prepared enzymes to evaluate convergence. Co-crystallized ligands were also docked in order to assess the docking performance through a re-docking strategy. Additionally, the specificity and sensitivity of the docking protocol, 15 compounds of diverse chemical nature, which have IC 50 < 4 µM, against COX-2 (3LN1) were compiled from ChEMBL [18]. Furthermore, 50 decoys per each active compound were compiled from the directory of useful decoys-enhanced (DUD-E) [19]. Thus, decoys and bioactive compounds were processed using the same docking protocol. The resulting data was then assessed through a receiver operating characteristic (ROC) and scores enrichment using the screening explorer webserver [20]. Once the docking parameters were validated, the docking simulations of test 2-arylbenzofurans were performed. Docking results were initially examined in PyMOL 2.3 and the resulting poses were top-ranked regarding the Vina scores (kcal/mol). Molegro Virtual Docker 6.0 (CLC Bio Company, Aarhus, Denmark) was secondly used to evaluate the best-poses of best-docked compounds as a rescoring strategy. RMSD of atomic positions of the best poses between replicates was calculated through an in-house bash script. Discovery Studio Client 16.1 (Biovia, San Diego, CA, USA) and Maestro 11.8 (Schrödinger, Cambridge, MA, USA) were used to visualize the respective 3D models and 2D residual interactions diagrams, respectively, using the pdbqt outputs for the highest-scored ligands.

Statistical Analysis
Descriptive and inferential statistics tests were carried out using R project software version 3.0.2 (R Foundation, Vienna, Austria). Principal component analysis (PCA) was also carried out in SIMCA 14.0 software (Umetrics, Umeå, Sweden) using the docking scores in order to observe plausible relationships with structural variations among test 2-arylbenzofurans.

Molecular Dynamics Simulations
The ligand's best poses from selected docked (ligands) and crystal (COX-2, 3LN1) structures were chosen as inputs for molecular dynamics (MD) simulations of ligand-enzyme complexes and the apoenzyme. Such simulations were carried out on the systems comprised of COX-2+celecoxib, COX-2 + 11, COX-2 + 29, COX-2 + 42, and COX-2 + 49 (obtained after molecular docking). These systems were individually subjected to an additional 150-ns conventional MD simulation for exploring the performance over the time of these ligands within the active site of the enzyme. These MD simulations were run in Gromacs 5.0.5 (open source, http://www.gromacs.org) on an Ubuntu 12.04 server, using NPT (constant pressure and temperature) and periodic boundary conditions, as previously reported [21]. Hence, docked ligands were prepared by adding hydrogen atoms in UCSF Chimera (UCSF, CA, USA) and the resulting pdb-file was uploaded to the atb server (http://compbio.biosci.uq.edu.au/atb/) to add the respective Gromo53a6 force field and get the itp-type topology file. COX-2 topologies were obtained in Gromacs using the Gromos53a6 force fields, due to the presence of the Heme group. The simple point charge (SPC) water model was then implemented for solvation in a triclinic box using a 1.0-nm margin distance. 0.10 M NaCl was added to the simulation systems and water molecules were randomly replaced until neutrality. An energy minimization through a 2000-steps steepest descent method was then used. NVT (constant volume and temperature) equilibration at 310 K during 500 ps, followed by NPT equilibration during 1000 ps using the Parrinello-Rahman method at 1 bar as a reference were done on the systems using position restraints. Lastly, solute position restraints were released and a production run for 90 ns was performed. Temperature and pressure were kept constant at 310 K and 1 bar, respectively. Coordinates were recorded in a 1 fs time step. Electrostatic forces were calculated using the particle-mesh Ewald (PME) method. Periodic boundary conditions were used in all simulations and covalent bond lengths were constrained by the LINear Constraint Solver (LINCS) algorithm.

Binding Free Energy Analyses
Binding-free energy was calculated using the g_mmpbsa tool [22] (open source drug discovery consortium (OSDD), New Dehli, India). This tool calculated components of the free energy of the protein-substrate binding (∆G Bind ) using the molecular mechanic/poisson-boltzmann surface area (MM/PBSA) method [23,24]. In this method, ∆G bind calculation between a protein and a ligand is carried out by ∆G bind = ∆H − T∆S ≈ ∆E MM + ∆G sol − T∆S, ∆E MM = ∆E internal + ∆E electrostatic − ∆E vdW , ∆G solv = ∆G elec solv + DG vdW solv , where the total gas phase energy on the binding of MM energy is ∆E MM , the free energy of solvation is ∆G solv , and the entropy contribution is T∆S. A Poisson-Boltzmann model was used to compute the electrostatic solvation energy in a continuum solvent. The derivation of a non-polar solvation energy term was computed as a solvent-accessible surface area (SASA). ∆E MM were calculated using the Lennard-Jones and Coulomb potential [23,24]. ∆G Bind was used to analyze the binding associations between COX-2 and selected ligands (i.e., celecoxib, 11, 29, 42, and 49) by decomposing the total binding-free energy into each residue. The binding energy calculations of the selected ligands were performed for 100 snapshots taken at an interval of 500 ps during the last stable 40-ns MD simulations.

COX-2 Inhibition Assay
The ability of selected compounds to inhibit bovine COX-2 was determined using an enzyme immunoassay (EIA) kit (Cayman Chemical Company, Ann Arbor, MI, USA), according to the reported method [8]. The activity was expressed as half-maximal inhibitory concentrations, IC 50 (in µM), obtained after non-linear regression, using the program GraphPad Prism version 5.00 (GraphPad Software, San Diego, CA, USA).

Comparative Molecular Field Analysis (CoMFA)
The best-docked poses of those test 2-arylbenzofurans experimentally evaluated against COX-2 were merged and aligned by means of particular tethers placed on benzofuran and aromatic rings moieties, using the molecular overlay tool included in the software Discovery Studio Client 16.1 (Biovia, San Diego, CA, USA). Compound 11 was selected as a template to select tethers. The aligned molecules set (in an sdf-file) was randomly divided into two subsets (training and test sets, corresponding to 70% and 30%, respectively). Comparative molecular field analysis (CoMFA) analysis was then performed using Open3DQSAR, using the standard protocol [25]. The experimental in-vitro COX-2 inhibition (expressed as IC 50 in M) was converted into a negative logarithmic form (pIC 50 ) and then used as an independent variable. The models were validated by leave-one-out (LOO) and leave-many-out (LMO) methods. The quality of the models after validation was evaluated on predicting the independent variable for the test set.

2-Arylbenzofuran-Related Chemical Space
The selection of test 2-arylbenzofuran-type compounds was based on three parameters: (1) natural origin, (2) enlarging the chemical space of an in-house collection (n = 26) of 2-arylbenzofurans, and (3) having an anti-inflammatory effect (preferably, inhibition of enzymatic activity or expression of COX-2). Hence, a targeted, custom-made library, comprising 58 2-arylbenzofurans, was then compiled. Common names of these compounds (1-58) are found in Table A1 and their structures are presented in Figure A2. The 2-arylbenzofuran-related chemical space was examined to characterize it and facilitate the data analysis according to structural fragments. Therefore, a structural similarity analysis between test compounds was performed using the FragFp descriptor available in DataWarrior software [17]. Compound 11 was selected as a reference to compute the similarity index due to its good COX-2 inhibitory activity, as previously reported [8]. The resulting similarity chart is presented in Figure 1.
After such an analysis, some clusters were shown, which involved similarity FragFp indexes between 0.93 and 0.28, according to the heatmap based on the structure similarity index between 0 (red) and 1 (light blue) ( Figure 1). Clusters were related to different classes of substituted benzofurans, which were then subdivided into moracin-type (including isopentyl-substituted and dihydrooxepine-containing) ( analysis clustered quantitatively the test compounds into six classes, which indicated these compounds have particular moieties/fragments to have further analysis. having an anti-inflammatory effect (preferably, inhibition of enzymatic activity or expression of COX-2). Hence, a targeted, custom-made library, comprising 58 2-arylbenzofurans, was then compiled. Common names of these compounds (1-58) are found in Table A1 and their structures are presented in Figure A2. The 2-arylbenzofuran-related chemical space was examined to characterize it and facilitate the data analysis according to structural fragments. Therefore, a structural similarity analysis between test compounds was performed using the FragFp descriptor available in DataWarrior software [17]. Compound 11 was selected as a reference to compute the similarity index due to its good COX-2 inhibitory activity, as previously reported [8]. The resulting similarity chart is presented in Figure 1. Figure 1. Similarity chart of 2-arylbenzofurans 1-58 after structure similarity analysis using the substructure fragment dictionary-based binary fingerprint descriptor (FragFp) ( [17]. Colors according to the heatmap based on the structure similarity index between 0 (red) and 1 (light blue). Clusters according Figure 1. Similarity chart of 2-arylbenzofurans 1-58 after structure similarity analysis using the substructure fragment dictionary-based binary fingerprint descriptor (FragFp) ( [17]. Colors according to the heatmap based on the structure similarity index between 0 (red) and 1 (light blue). Clusters according to: (a) moracin-type, (b) 2H-pyran-and benzofuro[6.5-b]furan-containing, (c) propyl(en)-susbtituted, (d) dihydropyran-containing, (e) furocyclohexadienone-type, and (f) licarin-type.

Molecular Docking Simulations
Molecular docking was chosen as a first-line discrimination of bioactives from the above-mentioned chemical space related to 2-arylbenzofuran compounds due to the capacity to simulate the binding of small compounds within the active site of enzymes, as a structure-based discrimination procedure for the in-silico prediction of putative competitors [26][27][28]. In this regard, some concerns should be considered using accurate parameters for suitable selection of sampling and scoring procedures during structure-based screening [29,30]. Therefore, the molecular docking of the co-crystallized ligands for each enzyme structure (ES1-ES9, Table 1) indicated a good performance of the present docking protocol since these re-docking calculations resulted into low conformational RMSD values (<0.5 Å), such as the docked/co-crystallized superposition of celecoxib presented in Figure 2a. Furthermore, the docking performance was also evaluated through a benchmarking approach from docking results of 15 known selective COX-2 inhibitors in comparison to the binding scores of a group of decoys (n = 750). Calculation of the area under the curve (AUC) from ROC curves was used to estimate the sensitivity and specificity of the docking protocol [31]. AUC and the Boltzmann-enhanced discrimination of the receiver operating characteristic (BEDROC) resulted in 0.908 and 0.827, respectively. Therefore, the validation of the docking protocol was considered successful (Figure 2b). This evaluation reached true positive identification over 92% of the active compounds involving a recognition within 10% of the test compounds ( Figure 2c).
Calculation of the area under the curve (AUC) from ROC curves was used to estimate the sensitivity and specificity of the docking protocol [31]. AUC and the Boltzmann-enhanced discrimination of the receiver operating characteristic (BEDROC) resulted in 0.908 and 0.827, respectively. Therefore, the validation of the docking protocol was considered successful (Figure 2b). This evaluation reached true positive identification over 92% of the active compounds involving a recognition within 10% of the test compounds (Figure 2c). Once the parameters were adequately validated, each geometrically-optimized structure 1-58 was docked into the active site of each COX-2 structures (ES1-ES9, Table 1) using AutoDock/Vina. These enzyme structures ES1-ES9 were retrieved from the PDB and their respective data are shown in Table 1. These M. musculus-derived COX-2 structures were selected to compare the performance of selective (S) and non-selective (NS) COX-2 inhibitors and 2-arylbenzofurans, as a discriminating Once the parameters were adequately validated, each geometrically-optimized structure 1-58 was docked into the active site of each COX-2 structures (ES1-ES9, Table 1) using AutoDock/Vina. These enzyme structures ES1-ES9 were retrieved from the PDB and their respective data are shown in Table 1. These M. musculus-derived COX-2 structures were selected to compare the performance of selective (S) and non-selective (NS) COX-2 inhibitors and 2-arylbenzofurans, as a discriminating initiative using the re-docking strategy, to assess the predictive potential of the simulation protocol. It was, therefore, required to retrieve COX-2 structures bound to common anti-inflammatory drugs (i.e., S/NS COX-2 inhibitors) and the natural ligand. In addition, in spite of the importance to include human COX-2, they were not used in the present docking study. Other MmCOX-2 PBD structures were discarded due to some issues related to sequence, residues, and a bound ligand. Such facts were taken as criterion to select the test COX-2 structures. Furthermore, the active site was adopted according to the information described in the referenced studies. After a detailed scrutiny of such information, five common residues involved in the binding site were initially defined (i.e., Arg106, Val335, Ser339, Tyr341, and Ser530). They were then used in the present docking study as flexible residues as well as those located at 4.0 Å from the co-crystallized ligand. These structures used for docking consist of very similar COX-2 protein sequences in a different solution when the enzyme was co-crystallized with a ligand or unbounded. These structures ES1-ES9 presented small alterations on the tertiary structure depending on the crystallization protocol. In addition, the presence of the inhibitor bound to the active site can also induce changes in the 3D-conformation of the binding site and modify the accessibility of some residues for ligand interaction. These facts were examined on superimposing the 9 chain-A 3D-structures of COX-2 after sequence alignment ( Figure 3), which exhibited slight differences in a primary and tertiary structure (i.e., sequence similarity 99.2%-99.8% and RMSD 0.35-0.50 Å (Figure 3a), respectively). The differences between these COX-2 sequences were found to be related to particular residues such as Asn34, Ala33, His90, His208, Gln310, and Lys333 ( Figure 3b), but the active site remained identical. Therefore, the same ligands docked into the same binding site of different crystal structures of the same enzyme is very important to statistically support the docking results.

Performance and Trends of the Docking Results
Docking results were examined by descriptive statistics to explore the trends and performance of the docking protocol, on one hand, and visualize and rationally select the best-docked compounds to discriminate putative binders of COX-2 on the other. After 10 replicates per calculation of benzofuran-related compounds 1-58 with nine crystal structures of COX-2 (ES1-ES9), the docking results were expressed as Vina Scores (in kcal/mol) and the best pose into the active site of each enzyme structure was also achieved. The mean docking scores, the standard deviation (SD) among replicates, and the RMSD (in Å) of atomic positions between the best poses after 10 replicates for each test benzofuran docked with each COX-2 structure are listed in Table A2. In general, the SD values were found to be within the 0.02-0.55 range, which indicates good reproducibility for docking calculations. However, the RMSD values exhibit high variability (i.e., 0.049-1.582 Å range). However, the most RMSD values were below 0.8 Å, which demonstrated good convergence in the calculated complexes formation. As expected, arachidonic acid (AA) exhibited the highest variability due to the long chain-derived conformations (14 rotatable bonds), but planar (more rigid) structures exhibited the lowest SD and RMSD values. 2-arylbenzofurans possessing an additional ring (e.g., pyran, furan, or oxepine) exhibited overall stronger docking scores (<−10 kj/mol). The previously mentioned facts were confirmed with the global mean Vina scores (GMVS) for the docking results with all COX-2 structures ( Table 2).
criterion to select the test COX-2 structures. Furthermore, the active site was adopted according to the information described in the referenced studies. After a detailed scrutiny of such information, five common residues involved in the binding site were initially defined (i.e., Arg106, Val335, Ser339, Tyr341, and Ser530). They were then used in the present docking study as flexible residues as well as those located at 4.0 Å from the co-crystallized ligand. These structures used for docking consist of very similar COX-2 protein sequences in a different solution when the enzyme was co-crystallized with a ligand or unbounded. These structures ES1-ES9 presented small alterations on the tertiary structure depending on the crystallization protocol. In addition, the presence of the inhibitor bound to the active site can also induce changes in the 3D-conformation of the binding site and modify the accessibility of some residues for ligand interaction. These facts were examined on superimposing the 9 chain-A 3D-structures of COX-2 after sequence alignment (Figure 3), which exhibited slight differences in a primary and tertiary structure (i.e., sequence similarity 99.2%-99.8% and RMSD 0.35-0.50 Å (Figure 3a), respectively). The differences between these COX-2 sequences were found to be related to particular residues such as Asn34, Ala33, His90, His208, Gln310, and Lys333 (Figure 3b), but the active site remained identical. Therefore, the same ligands docked into the same binding site of different crystal structures of the same enzyme is very important to statistically support the docking results.  Additionally, docking scores were found to be adequately distributed and reproducible. The relative standard deviation (RSD) percentages resulted in the 2.7%-33.2% range, but only 16.3% of the RSD values were above 10% and 3.6% of the RSD values above 20%. These facts are a good indication of the reproducibility. It is well-recognized that the docking procedure has a limitation related to distinguishing the correct pose among the generated scoring positions. However, the difference in the best poses among replicates were considered to be convergent since 76.3% of the RMSD values were below 0.6 Å. A similar trend was observed with inhibitors and the respective co-crystallized structure after re-docking calculations (RMSD < 0.5 Å). Table 2. Global mean values of Vina scores, relative standard deviation percentage (%RSD), and root mean square deviation (RMSD) (in Å) of 2-arylbenzofurans docked with cyclooxygenase-2 (COX-2) structures. Regarding the behavior of Vina scores among COX-2 structures (ES1-ES9), a similar distribution was found as shown in the respective boxplot (Figure 4a), which was built from the resulting mean docking scores. Mean values per COX-2 structure were similar for ES1-ES5, as a first group, and ES6, ES8, and ES9 as a second one. ES7 exhibited the lowest mean Vina score among enzymes and the lowest docking score (−12.7 ± 0.2 kcal/mol) corresponding to celecoxib (cel). These facts mean that the particular behavior of the ligands into each binding pocket was slightly different, which can be understood as a consequence of structural adaptations of the active site by the presence of the co-crystallized ligands. This issue should be considered a very important factor in the moment to get more reliable results during docking studies.

Performance and Trends of the Docking Results
Biomolecules 2020, 10, 176 9 of 23 ES8, and ES9 as a second one. ES7 exhibited the lowest mean Vina score among enzymes and the lowest docking score (−12.7 ± 0.2 kcal/mol) corresponding to celecoxib (cel). These facts mean that the particular behavior of the ligands into each binding pocket was slightly different, which can be understood as a consequence of structural adaptations of the active site by the presence of the cocrystallized ligands. This issue should be considered a very important factor in the moment to get more reliable results during docking studies. Trends within docking results were also inspected by multivariate statistics in order to identify the plausible relationship between docking scores and structural variations of 2-arylbenzofurans. A PCA was performed using the Vina scores dataset. The first two components explained the 82.3% of the total variance with good model fitting (R 2 = 0.823 and Q 2 = 0.778 for the principal component 1 (PC1) and PC2). The resulting score plot (PC1 vs. PC2, Figure 4b) demonstrated good differentiation of the compounds by variance according to the Vina scores along COX-2 structures. Hence, this unsupervised analysis clustered the compounds in four groups, which is shown by discriminating across PC1 depending on mean Vina scores (66.2% explained variance). Group 1 (green) consisted of the best-docking 2-arylbenzofurans (such as 11, 28, 29, 30, 40, 42, 49, 51, and 54) with a comparable profile to the selective COX-2 inhibitors celecoxib (cel) and SC558. 2-arylbenzofurans possessing an additional cycle (such as pyran, furan, and oxepine) were the predominant structural feature in group 1. In contrast, group 2 (dark blue) clustered the weakest-docking 2-arylbenzofurans with profiles related to those of arachidonic acid (AA), indomethacin, and naproxen. This group was characterized by unsubstituted moracin-type compounds as well as some dihydrobenzofurans. Groups 3-4 (yellow and red) exhibited their separation by dual incidence of PC1 and PC2, whose separation along PC2 was mainly dependent by the scoring behavior with ES6 and ES8. These results indicated that the PCA led to discriminated bioactive compounds, according to the molecular docking performance.

Selection of the Best-Docked Compounds
An examination of the previously mentioned docking results for compounds 1-58 (Table 2) was performed to select the best-docked compounds as putative COX-2 binders. Hence, from the global mean Vina scores (GMVS < −10.5 kj/mol)), best-pose convergence (RMSD < 0.5) and scoring reproducibility (%RSD <6%) along nine COX-2 structures, the compounds 9′-nor-7, 8-dehydro-isolicarin-B (11), moracin H (29) psoralidin (42), and maximol (49) were established as the best-docked compounds among the test set. These selected compounds were, therefore, clustered in the green group of the PCA-derived score plot (Figure 4b), which exhibits the lowest-scoring behavior. Other compounds in this group were not selected Trends within docking results were also inspected by multivariate statistics in order to identify the plausible relationship between docking scores and structural variations of 2-arylbenzofurans. A PCA was performed using the Vina scores dataset. The first two components explained the 82.3% of the total variance with good model fitting (R 2 = 0.823 and Q 2 = 0.778 for the principal component 1 (PC1) and PC2). The resulting score plot (PC1 vs. PC2, Figure 4b) demonstrated good differentiation of the compounds by variance according to the Vina scores along COX-2 structures. Hence, this unsupervised analysis clustered the compounds in four groups, which is shown by discriminating across PC1 depending on mean Vina scores (66.2% explained variance). Group 1 (green) consisted of the best-docking 2-arylbenzofurans (such as 11, 28, 29, 30, 40, 42, 49, 51, and 54) with a comparable profile to the selective COX-2 inhibitors celecoxib (cel) and SC558. 2-arylbenzofurans possessing an additional cycle (such as pyran, furan, and oxepine) were the predominant structural feature in group 1. In contrast, group 2 (dark blue) clustered the weakest-docking 2-arylbenzofurans with profiles related to those of arachidonic acid (AA), indomethacin, and naproxen. This group was characterized by unsubstituted moracin-type compounds as well as some dihydrobenzofurans. Groups 3-4 (yellow and red) exhibited their separation by dual incidence of PC1 and PC2, whose separation along PC2 was mainly dependent by the scoring behavior with ES6 and ES8. These results indicated that the PCA led to discriminated bioactive compounds, according to the molecular docking performance.

Selection of the Best-Docked Compounds
An examination of the previously mentioned docking results for compounds 1-58 (Table 2) was performed to select the best-docked compounds as putative COX-2 binders. Hence, from the global mean Vina scores (GMVS < −10.5 kj/mol)), best-pose convergence (RMSD < 0.5) and scoring reproducibility (%RSD <6%) along nine COX-2 structures, the compounds 9 -nor-7, 8-dehydro-isolicarin-B (11), moracin H (29) psoralidin (42), and maximol (49) were established as the best-docked compounds among the test set. These selected compounds were, therefore, clustered in the green group of the PCA-derived score plot (Figure 4b), which exhibits the lowest-scoring behavior. Other compounds in this group were not selected because they did not meet completely the previously mentioned features (i.e., GMVS, RMSD, and %RSD). In addition, Molegro Virtual Docker (MVD) was also used to simulate the intermolecular interaction as a rescoring strategy. An identical trend was obtained, so the best-docked compounds used in Vina exhibited the best MolDock scores (i.e., cel = −163.2, 49 = −156.9, 42 = −156.6, 11 = −139.1, 29 = −136.2), which is another indication of the good performance of the docking protocol and the filtering by descriptive and multivariate statistics. These four compounds were used for the first time to simulate their interaction within the active site of COX-2 using a validated docking protocol. However, compound 11 was previously evaluated against COX-2 enzyme (affording an IC 50 of 3.32 µM) [8] and compound 42 was reported to exhibit good ability to reduce the COX-2 expression at 50 µM [32].

Binding Mode and Residual Interactions of the Best-Docked Compounds
Binding modes of the best-docked 2-arylbenzofurans with the COX-2 structure ES7 (PDB code: 3LN1) was examined looking for important/crucial interactions between active site residues and ligand moieties, according to the individual 2D-residual interaction diagrams ( Figure 5).

Binding Mode and Residual Interactions of the Best-Docked Compounds
Binding modes of the best-docked 2-arylbenzofurans with the COX-2 structure ES7 (PDB code: 3LN1) was examined looking for important/crucial interactions between active site residues and ligand moieties, according to the individual 2D-residual interaction diagrams ( Figure 5). Analysis of these diagrams was useful to delineate some insights into the binding mode of the simulated complexes, which showed some important hydrophobic and polar interactions. These interactions were taken as key contacts and an indication of the importance of the presence of the hydroxiaryl moiety for interacting with the respective active site residues to stabilize the simulated enzyme-ligand complexes. Particularly, celecoxib exhibited H-bonds between Phe504, Arg499 (as donors), Analysis of these diagrams was useful to delineate some insights into the binding mode of the simulated complexes, which showed some important hydrophobic and polar interactions. These interactions were taken as key contacts and an indication of the importance of the presence of the hydroxiaryl moiety for interacting with the respective active site residues to stabilize the simulated enzyme-ligand complexes. Particularly, celecoxib exhibited H-bonds between Phe504, Arg499 (as donors), and Gln178 (as acceptor) and sulfonamide as well as Arg106 (as donor) with 1H-pyrazol-2-ium moiety. In contrast, other chemical interactions were found on analyzing the binding modes of the lowest-scored poses for the four strongest-docking compounds (i.e., 11, 29, 42, and 49; Figure 5). These compounds resulted in the lowest mean docking scores (−10.48 to −11.26 kcal/mol range) and were found to be well-positioned into the COX-2 active site. However, they involve different orientations depending on the type of interactions. In this regard, residues Tyr341, Gln178, and Ser339 were found to be key polar contacts to stabilize the enzyme ligand complex for those compounds possessing a resorcinyl substituent, through a π-π-stacking for 29 and acting as an H-donor for 49. Furthermore, the p-hydroxyphenyl group in 49 interacted with Trp373 via π-π-stacking. Psoralidin (42) exhibited another binding mode involving His75 as an H-donor to its phenolic OH at benzofuran moiety. Ala513 also interacted with phenolic OH at aryl moiety in 42. In addition, the methylenedioxy moiety in 11 was important to orient the molecule toward Arg499, acting as an H-donor acceptor. The combination of these structural features of such benzofurans could serve as an important starting point to design a novel series of COX-2 inhibitors.

Stability and Chemical Reactivity of the Best-Docked 2-Arylbenzofurans
Frontier molecular orbitals HOMO and LUMO can reflect the successful biological interaction of ligands within the binding pocket of proteins [33]. In this sense, the DFT level-derived calculations of HOMO and LUMO energies were a suggestion of kinetic stability and chemical reactivity of these best docked compounds. The resulting B3LYP/6-311++G (2df, 2p) models are presented in Figure 6. These compounds exhibited lower HOMO-LUMO energy gaps to that of celecoxib and other 2-arylbenzofurans, which infer high stability but also high reactivity. Pearson's correlation between COX-2 docking scores and HOMO-LUMO energy gaps for 2-arylbenzofurans 1-58 was also evaluated. A positive correlation among scores and energy gaps was found (0.812, p < 0.001), which indicated that the self-reactivity of compounds might influence the binding mode. Hence, among the best-docked compounds, 11 and 29 exhibited higher Vina scores and HOMO-LUMO energy gaps, whereas 42 and 49 showed an opposite behavior ( Figure 6). This relationship can be rationalized since these energy gaps and eigenvalues defined the biological activity of ligands via high reactivity and good stability. Therefore, enzyme inhibition by these best-docked 2-arylbenzofurans might be afforded through an electron donating ability of test compounds with the enzyme-binding site pocket of COX-2. and Gln178 (as acceptor) and sulfonamide as well as Arg106 (as donor) with 1H-pyrazol-2-ium moiety. In contrast, other chemical interactions were found on analyzing the binding modes of the lowest-scored poses for the four strongest-docking compounds (i.e., 11, 29, 42, and 49; Figure 5). These compounds resulted in the lowest mean docking scores (−10.48 to −11.26 kcal/mol range) and were found to be well-positioned into the COX-2 active site. However, they involve different orientations depending on the type of interactions. In this regard, residues Tyr341, Gln178, and Ser339 were found to be key polar contacts to stabilize the enzyme ligand complex for those compounds possessing a resorcinyl substituent, through a π-π-stacking for 29 and acting as an H-donor for 49. Furthermore, the phydroxyphenyl group in 49 interacted with Trp373 via π-π-stacking. Psoralidin (42) exhibited another binding mode involving His75 as an H-donor to its phenolic OH at benzofuran moiety. Ala513 also interacted with phenolic OH at aryl moiety in 42. In addition, the methylenedioxy moiety in 11 was important to orient the molecule toward Arg499, acting as an H-donor acceptor. The combination of these structural features of such benzofurans could serve as an important starting point to design a novel series of COX-2 inhibitors.

Stability and Chemical Reactivity of the Best-Docked 2-Arylbenzofurans
Frontier molecular orbitals HOMO and LUMO can reflect the successful biological interaction of ligands within the binding pocket of proteins [33]. In this sense, the DFT level-derived calculations of HOMO and LUMO energies were a suggestion of kinetic stability and chemical reactivity of these best docked compounds. The resulting B3LYP/6-311++G (2df, 2p) models are presented in Figure 6. These compounds exhibited lower HOMO-LUMO energy gaps to that of celecoxib and other 2-arylbenzofurans, which infer high stability but also high reactivity. Pearson's correlation between COX-2 docking scores and HOMO-LUMO energy gaps for 2-arylbenzofurans 1-58 was also evaluated. A positive correlation among scores and energy gaps was found (0.812, p < 0.001), which indicated that the self-reactivity of compounds might influence the binding mode. Hence, among the best-docked compounds, 11 and 29 exhibited higher Vina scores and HOMO-LUMO energy gaps, whereas 42 and 49 showed an opposite behavior ( Figure 6). This relationship can be rationalized since these energy gaps and eigenvalues defined the biological activity of ligands via high reactivity and good stability. Therefore, enzyme inhibition by these best-docked 2-arylbenzofurans might be afforded through an electron donating ability of test compounds with the enzyme-binding site pocket of COX-2.

Molecular Dynamics Studies of the Best-Docked 2-Arylbenzofurans
In order to extend the information on the binding mode of the best docked compounds, 150-ns MD simulations were separately accomplished using the COX-2 (ES7) alone (apoenzyme) and docked distinctly with compounds 11, 29, 42, 49, and celecoxib. Ligand-enzyme trajectories for resulting

Molecular Dynamics Studies of the Best-Docked 2-Arylbenzofurans
In order to extend the information on the binding mode of the best docked compounds, 150-ns MD simulations were separately accomplished using the COX-2 (ES7) alone (apoenzyme) and docked distinctly with compounds 11, 29, 42, 49, and celecoxib. Ligand-enzyme trajectories for resulting complexes were examined through the variation of geometric properties over the time. Thus, RMSD of the COX-2 backbone reflected the receptor frame stability by computing the time-dependent distance (Å) among different positions of the atom set (Figure 7a). Apoenzyme exhibited a normal evolution during the simulation but revealed a slight perturbation at 30 and 50 ns (RMSD 0.50-0.65 nm). A normal stabilization was then reached throughout the rest of the MD simulation (RMSD 0.63 nm). The COX-2·celecoxib complex evolved typically but showed a late perturbation at 80 ns (until RMSD 0.99 nm). Although compounds 11 and 49 exhibited the most-perturbed profile, both complexes (with celecoxib and 11) attained stability after 90 ns and 49 after 110 ns. The steady progress was maintained over the remaining time for these three ligands. On the other hand, compounds 29 and 42 exhibited the least-perturbed MD performance and the evolution was found to be similar to that of apoenzyme, which indicates good geometric properties on interacting with COX-2. The computed variations in RMSD values for these complexes showed that the COX-2 structure is affected differentially by the interaction with each ligand, which achieved a stable condition at the end of each MD simulation.  Fluctuations of the Cα atomic positions for each residue of the enzyme was explored through the root mean square fluctuation (RMSF). It was used to scrutinize the flexibility and secondary structure of the COX-2 enzyme under binding with the best-docked 2-arylbenzofurans. Thus, lower RMSF values implied constrained regions whereas higher RMSF values denoted more flexibility. All Fluctuations of the Cα atomic positions for each residue of the enzyme was explored through the root mean square fluctuation (RMSF). It was used to scrutinize the flexibility and secondary structure of the COX-2 enzyme under binding with the best-docked 2-arylbenzofurans. Thus, lower RMSF values implied constrained regions whereas higher RMSF values denoted more flexibility. All MD-simulated complexes showed similar performance (Figure 7b) involving fluctuations within the range of 0.2 and 1.6 nm. Such alterations were found to be near the formerly recognized crucial interactions (e.g., Arg106, Phe186, Leu292, Val330, Leu377, and Ser516). However, in comparison to the apoCOX-2, the binding with test ligands did not substantially modify the flexibility of these regions to maintain an overall steady state (RMSF differences > 0.2 nm), excepting the most-fluctuating residues around Arg106. Therefore, the COX-2 inhibition mode of compounds 11, 29, 42, and 49 might be interpreted through a complex stabilization, which is in agreement with the previously reported MD study for other known COX-2 inhibitors [34].

Binding-Free Energy Calculations of the Best-Docked 2-Arylbenzofurans
The binding-free energies (∆G bind ) for the best-docked compounds (i.e., 11, 29, 42, 49, and celecoxib) during the interaction of COX-2 for the last 40 ns of MD trajectory was estimated by the MM/PBSA approach to evaluate the global stability of the resulting ligand-enzyme complexes. The calculated binding energies are presented in Table 3. All five ligands exhibited negative binding energies. However, 42 showed comparable binding energy to that of celecoxib (−195.7 ± 4.5 kJ/mol vs. −197.8 ± 4.5 kJ/mol, respectively) and significant differences to that of other ligands, which justify the attained docking performance. The main contribution to the binding energy was due to van der Waal (∆E vdW ) energies (< −205 kJ/mol), which resulted in similar values among test compounds, while the contribution of electrostatic energy exhibited higher differences between them. Celecoxib and 49 exhibited the lowest electrostatic energy whereas 11 exhibited the highest one. In addition, the contribution of the polar solvation energy was found to be more unfavorable for 49 and celecoxib. The solvent accessible surface area (SASA) energy was similar for test ligands, even though compound 49 exhibited the highest ∆G sasa . Therefore, according to these results, electrostatic and polar solvation energy contributions could rationalize the difference in the binding mode of test ligands, which plays hydrophobic interactions as an important role in stabilization and even binding of test 2-arylbenzofurans within the active site of COX-2. As a consequence, non-polar electrostatic interactions could be implied as the main driving force for the molecular recognition of COX-2 by test ligands. This fact was then confirmed after per-residue decomposition of the total binding energy of the simulated COX-2 ligand complexes. In this context, most binding energy contributing residues were found to be different depending on the ligand, as presented in Table 4. Regarding ligands, 42 exhibited similar binding energy to that of celecoxib (<−80 kJ/mol), which has the highest values among test ligands. Regarding residues, Arg106, Leu338, and Val509 were found to be those residues that contribute the most (<−8 kJ/mol) to the binding energy for the COX-2·celecoxib complex (as described previously [34]), whereas Arg106 and Arg499 were observed for COX-2·11, Val335 for COX-2·29, Val335, Val509, and Ala513 for COX-2·42, and Glu510 and Ala513 for the COX-2·49 complex. Ser339 (<−2 kJ/mol) and Val509 (<−6 kJ/mol) were the common residues that contributed to the binding energy for all test ligands. Such residues with higher binding energy contributions represented those contact points for electrostatic interactions, as described in molecular docking and Cα atom fluctuations. Similar information was described for the binding mode of stilbene analogs within the active site of COX-2 [35].

In-Vitro COX-2 Activity of an in-House Collection of 2-Arylbenzofurans
The information on the potential as COX-2 inhibitors for the best-docked compounds (i.e., 11, 29, 42, 49) was expanded by measuring the IC 50 (in µM) against COX-2 of an in-house library of 26 2-arylbenzofurans, which several of them were evaluated in a previous study [8]. The resulting IC 50 values are presented in Table 5. Test compounds exhibited IC 50 values in the range of 438 and 0.752 µM. In general, furocyclohexadienone and dihydrobenzofuran-containing compounds exhibited the highest IC 50 values, except compound 49, whose IC 50 was 1.25 µM, which can be attributed to the resorcinyl and styryl substitutions at C3 and C5, respectively. These moieties influenced substantially to the binding of 49 since some interactions were favored with hydrophobic (i.e., styryl) and polar (i.e., resorcinyl) regions of the active site of COX-2. No significant differences were found for the IC 50 values of compounds 11 and 29, whose performances in molecular docking and molecular dynamics were very similar. The most active compound among the compound set was 42 (IC 50 = 0.752 µM), which exhibited an interesting docking behavior (mean Vina score −11.25 ± 0.59 kcal/mol), the highest stability and reactivity (∆E HOMO-LUMO = 3.77 eV), and an excellent performance over the time for stabilizing the COX-2·42 complex. This predicted compound was previously reported to exhibit anti-inflammatory properties as a dual inhibitor of the expression of COX-2 and 5-LOX [32]. These results can be understood as an adequate validation of the molecular docking-derived discrimination to find COX-2 inhibitors based on 2-arylbenzofurans. Previously mentioned most-active compounds also exhibited favorable ADMET properties (calculated using the FAF-Drugs4 web server [36]), which is afforded as accepted since the overall result with no Lipinski's rule violation, good solubility, and oral bioavailability (Table A3).

Comparative Molecular Field Analysis (CoMFA)
A CoMFA was applied to correlate molecular interaction fields (MIF) and COX-2 inhibition in order to consistently rationalize the inhibitory activity of 2-arylbenzofurans regarding electrostatic and steric properties [38] as well as describing and predicting some structural requirements for COX-2 inhibition by 2-arylbenzofurans. For this purpose, the best-docked pose of compound 11 was used as a template to outline tethers and overlay the test compounds (Table 5), due to the structural features. Aligned structures and experimental data of COX-2 inhibition (as pIC 50 values) were segmented into a training set (70%) to create the respective CoMFA model and a test set (30%) for the external validation [38]. MIFs were computed by means of electrostatic and steric probes. Partial least squares (PLS) regression was performed (employing up to five PLS components). Linear relationships between MIF fluctuations according to variations of the experimental pIC 50 were accomplished. Thus, the best model achieved good correlation between MIF values and experimental pIC 50 of test compounds using four PLS components, which comprises a correlation coefficient r 2 = 0.816, a F-test = 162.321, cross-validated LOO coefficient q 2 LOO = 0.766, and a cross-validated LMO coefficient q 2 LMO = 0.701. These parameters were found to be enough for a statistically robust, predictive model [39]. The adequate correlation was examined through a Y-scrambling procedure [40]. Hence, 30 scramblings and 10 runs were employed. After that, no correlation appeared since the models dropped significantly (R scr 2 and Q scr 2 < 0.3). This fact suggested the model was not reached as a consequence of a coincidental correlation. Therefore, this CoMFA-based model predicted the activity for the entire compound set, expressed as pIC 50 (pred) ( Table 5), which demonstrated a reasonable correlation for both training and test datasets. Steric and electrostatic field outputs (stdev * coeff) reached values of 31.1% and 37.4%, respectively. Figure 8 showed the corresponding translated contour surfaces from the field contributions. The electrostatic field contour map (Figure 8a) displayed the positively and negatively-charged regions to favor COX-2 inhibition as blue and red contours, respectively. Accordingly, phenolic hydroxyl groups at C4 and C7, resorcinyl moiety, and electron-withdrawing groups at C6 enhanced the COX-2 activity to interact with some residues located in the polar region of the active site of this enzyme (e.g., Hiz75, Arg499, and Gln178), as observed for compounds 29, 42, and 49. Alkyl (such as allyl, styryl, and prenyl) and methoxyl substitutions on 2-aryl moiety also favored the COX-2 inhibition (e.g., compounds 11 and 42) to interact positively or negatively with residues in the hydrophobic zone of the COX-2 active site (i.e., Val102, Leu103, and Leu345). On the other side, the enzymatic inhibitory activity might also be enhanced/decreased by the presence of bulky groups at benzofuran and 2-aryl moieties, as observed in the steric field contour maps (Figure 8b). Unfavorable and favorable zones by steric effects were then depicted as yellow and green contours, respectively. COX-2 inhibitory activity could be improved if a bulky steric group is located on C3′, C5, and C7, such as substituted prenyl groups and dihydropyran ring, as observed in most active compounds 29 and 42. In contrast, a bulky group in C2′ (e.g., hydroxypropenyl) and C3 (e.g., phenyl group) negatively influenced the enzymatic inhibition. Consequently, the CoMFA model indicated that the particular presence of both bulky and electronegative substituents (oriented sterically and electrostatically toward polar and hydrophobic zones within the active site) could be considered as important structural features to improve the COX-2 inhibitory activity by 2-arylbenzofurans.

Conclusions
An in-silico study was primarily performed on a custom-made library comprising 58 benzofurancontaining compounds of natural origin in order to discriminate binders and non-binders of COX-2. Thus, those 2-arylbenzofurans, having an additional ring (such as oxepine) or bulky substitutions in benzofuran or aryl moieties, can be considered as potential hits for COX-2 inhibition due to the favorable characteristics for an adequate binding within the active site of this enzyme. In this regard, the best-docked compounds exhibited the best performance in MD simulations, the most-favorable binding energies, and the lowest in-vitro IC50 values for COX-2 inhibition. Consequently, the molecular docking behavior of 2-arylbenzofurans within the active site of this enzyme might be used as an appropriate strategy to predict/analyze their binding mode and potentiality as COX-2 inhibitors. Lastly, the findings through the CoMFA model demonstrated a structure-activity relationship, which shows some crucial structural requirements for enhanced enzyme inhibition. These structural features must be oriented sterically and electrostatically toward both polar and hydrophobic zones within the active site to promote crucial interactions. These oriented interactions are the main driving force for the molecular recognition of COX-2 by test ligands. In this context, four hits (i.e., 11, 29, 42, and 49) are then proposed as lead structures for the development of COX-2 inhibitors. Therefore, such information should be conserved during further studies by focusing on the development of COX-2 inhibitors based on 2-arylbenzofurans.  On the other side, the enzymatic inhibitory activity might also be enhanced/decreased by the presence of bulky groups at benzofuran and 2-aryl moieties, as observed in the steric field contour maps (Figure 8b). Unfavorable and favorable zones by steric effects were then depicted as yellow and green contours, respectively. COX-2 inhibitory activity could be improved if a bulky steric group is located on C3 , C5, and C7, such as substituted prenyl groups and dihydropyran ring, as observed in most active compounds 29 and 42. In contrast, a bulky group in C2 (e.g., hydroxypropenyl) and C3 (e.g., phenyl group) negatively influenced the enzymatic inhibition. Consequently, the CoMFA model indicated that the particular presence of both bulky and electronegative substituents (oriented sterically and electrostatically toward polar and hydrophobic zones within the active site) could be considered as important structural features to improve the COX-2 inhibitory activity by 2-arylbenzofurans.

Conclusions
An in-silico study was primarily performed on a custom-made library comprising 58 benzofuran-containing compounds of natural origin in order to discriminate binders and non-binders of COX-2. Thus, those 2-arylbenzofurans, having an additional ring (such as oxepine) or bulky substitutions in benzofuran or aryl moieties, can be considered as potential hits for COX-2 inhibition due to the favorable characteristics for an adequate binding within the active site of this enzyme. In this regard, the best-docked compounds exhibited the best performance in MD simulations, the most-favorable binding energies, and the lowest in-vitro IC 50 values for COX-2 inhibition. Consequently, the molecular docking behavior of 2-arylbenzofurans within the active site of this enzyme might be used as an appropriate strategy to predict/analyze their binding mode and potentiality as COX-2 inhibitors. Lastly, the findings through the CoMFA model demonstrated a structure-activity relationship, which shows some crucial structural requirements for enhanced enzyme inhibition. These structural features must be oriented sterically and electrostatically toward both polar and hydrophobic zones within the active site to promote crucial interactions. These oriented interactions are the main driving force for the molecular recognition of COX-2 by test ligands. In this context, four hits (i.e., 11, 29, 42, and 49) are then proposed as lead structures for the development of COX-2 inhibitors. Therefore, such information should be conserved during further studies by focusing on the development of COX-2 inhibitors based on 2-arylbenzofurans.
Funding: This study was financially supported by The Vicerrectoría de Investigaciones at Universidad Militar Nueva Granada -Validity 2018, through the Project IMP-CIAS-2924 financially supported this study.

Conflicts of Interest:
The author declares no conflict of interest. The funders had no role in the design of the study, in the collection, analyses, or interpretation of data, in the writing of the manuscript, or in the decision to publish the manuscript.
Appendix A Table A1. Common names of test 2-arylbenzofuran-related compounds.

No.
Common Name Reference No.        Table A3. Absorption, distribution, metabolism, and excretion-toxicity (ADMET) properties of best-docked compounds within the active site of COX-2.